Info

Diese Frage ist geschlossen. Öffnen Sie sie erneut, um sie zu bearbeiten oder zu beantworten.

Problem with Finite Difference Scheme

1 Ansicht (letzte 30 Tage)
Danila Zharenkov
Danila Zharenkov am 6 Nov. 2013
Geschlossen: MATLAB Answer Bot am 20 Aug. 2021
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

Antworten (0)

Diese Frage ist geschlossen.

Produkte

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by