Problem with Finite Difference Scheme

Danila Zharenkov

2013 11 6

0 个回答

更新时间：2021 8 20

4 次查看（30 天）

Follow Question

信息

此问题已关闭。请重新打开它进行编辑或回答。

显示更早的评论

在 MATLAB Online 中打开

0 个投票

Hello, I'm trying to build the finite difference scheme for this task

I build the finite difference scheme like this

The main problem is boundary conditions for theta, because boundary values on its edges don't depend of other counted values. Here is my code, maybe I made a mistake in it?

clc;
clear all;
%constants
r0=0.1;
a=1.01;
b=1;
mu=1;
c=100;
%step-size for r
n=100;
r_max=1;
hr=r_max/n;
r=0:hr:r_max;
nr=max(size(r));
%step-size for theta
m=200;
th_max=2*pi;
hth=th_max/m;
th=0:hth:th_max;
nth=max(size(th));
%step-size for t
l=100;
t_max=5;
ht=t_max/l;
time=0:ht:t_max;
nt=max(size(time));
u=zeros(nr,nth,nt);
v=zeros(nr,nth,nt);
%Initialization
for i=1:nr
    for j=1:nth
        if ((i*hr<=r0)&&(0<=th(j))&&(th(j)<=pi/6))
           u0=0; 
        elseif ((i*hr<=r0)&&(pi/6<th(j))&&(th(j)<2*pi))
           u0=1;
        elseif ((i*hr>r0)&&(0<=th(j))&&(th(j)<=2*pi))
           u0=0; 
        end
        u(i,j,1)=u0;
    end
end
for i=1:nr
    for j=1:nth
        if ((r(i)<=r0)&&(0<=th(j))&&(th(j)<=pi/6))
           v0=1; 
        elseif ((r(i)<=r0)&&(pi/6<th(j))&&(th(j)<2*pi))
           v0=0;
        elseif ((r(i)>r0)&&(0<=th(j))&&(th(j)<=2*pi))
           v0=0; 
        end
        v(i,j,1)=v0;
    end
end
for t=1:nt-1
      for i=2:nr-1
          for j=2:nth-1
          u(i,j,t+1)=u(i,j,t)+ht*((u(i+1,j,t)-2*u(i,j,t)+u(i-1,j,t))/hr^2+(1/(i*hr)))*((u(i+1,j,t)-u(i,j,t))/hr)+(1/(i*hr)^2)*((u(i,j+1,t)-2*u(i,j,t)+u(i,j-1,t)/hth)+u(i,j,t)*(1-u(i,j,t)-c*v(i,j,t)));
          v(i,j,t+1)=v(i,j,t)+ht*((v(i+1,j,t)-2*v(i,j,t)+v(i-1,j,t))/hr^2+(1/(i*hr)))*((v(i+1,j,t)-v(i,j,t))/hr)+(1/(i*hr)^2)*((v(i,j+1,t)-2*v(i,j,t)+v(i,j-1,t)/hth)+a*v(i,j,t)*(1-c*u(i,j,t)-b*v(i,j,t)));
          %Boundary
          u(1,j,t)=u(2,j,t);
          u(nr,j,t)=u(nr-1,j,t);
          v(2,j,t)=v(1,j,t);
          v(nr,j,t)=v(nr-1,j,t);
          u(i,nth,t)=u(i,1,t);
          v(i,nth,t)=v(i,1,t);
          end
      end
  end