16 views (last 30 days)

Show older comments

I tried several methods, but i couldn't find the solution. for instance, i used Crunk-Nicelson finite difference method like following script but i don't know how can i apply the secend eq.(ODE) inside the matrix. the big problem here is that each incerment is too small.

% The Crank Nicolson Method

clc, clear, close

eps=0.5; tor=3; lf=5*10^-6;%microM dm=60*8.6*10^-6; %m2/s kg= 1.87*60; % m/s area=0.6 ;% m^-1

cin=10; % ki=1*60; %mg/m3s Ki=0.24; %m3/mg miu=0.3*10^-6; %microM^-1 beta=0.5; iif=16.5*10;%mw/cm2

Kw=4.9*10^-4;% m^3/mg cw=1000; %mg/m3 ux=0.011*60; %m/s

% syms x % d=int((exp(-miu*lf*x)^beta)/lf,0,lf);

%%%%%%%%%%%%%%%%%%parameter

f=iif^beta; g=ki*Ki; h=Kw*cw; alfa=eps*dm/tor;

tf=180; nx=100; dx=lf/nx; nt=2000; dt=tf/nt;

% --- Constant Coefficients of the tridiagonal system

c = alfa/(2*dx^2)+ux/(4*dx); % Subdiagonal: coefficients of u(i-1)

b = alfa/(2*dx^2)-ux/(4*dx); % Super diagonal: coefficients of u(i+1)

a = 1/dt+b+c; % Main Diagonal: coefficients of u(i)

% Boundary conditions and Initial Conditions

Uo(1)=20; Uo(2:nx)=0;

Un(1)=20; Un(nx)=0;

% Store results for future use

UUU(1,:)=Uo;

% Loop over time

for k=2:nt

for ii=1:nx-2

if ii==1

d(ii)=c*Uo(ii)+(1/dt-b-c)*Uo(ii+1)+b*Uo(ii+2)+c*Un(1);%-f*g*Uo(ii+1)/(1+h+Ki*Uo(ii+1));

elseif ii==nx-2

d(ii)=c*Uo(ii)+(1/dt-b-c)*Uo(ii+1)+b*Uo(ii+2)+b*Un(nx);%-f*g*Uo(ii+1)/(1+h+Ki*Uo(ii+1));

else

d(ii)=c*Uo(ii)+(1/dt-b-c)*Uo(ii+1)+b*Uo(ii+2);%-f*g*Uo(ii+1)/(1+h+Ki*Uo(ii+1));

end

end % note that d is row vector

%%%%%%%%%%%%%%%%%

% Transform a, b, c constants in column vectors:

bb=b*ones(nx-3,1);

cc=c*ones(nx-3,1);

aa=a*ones(nx-2,1);

% Use column vectors to construct tridiagonal matrices

AA=diag(aa)+ diag(-bb,1)+ diag(-cc,-1);

% Find the solution for interior nodes i=2,3,4,5

% UU=AA\d';

UU=inv(AA)*d';

% Build the whole solution as row vector

Un=[Un(1),UU',Un(nx)];

Un(nx)=Un(nx-1);

UUU(k,:)=Un;

Uo=Un;

end

t=[0:dt:tf-dt];

uf=UUU(:,end);

plot(t,uf)

Torsten
on 13 Feb 2018

Take a look at the answer provided here:

https://de.mathworks.com/matlabcentral/answers/371313-error-in-solving-system-of-two-reaction-diffusion-equations

The problem is very similar to yours.

Best wishes

Torsten.

Torsten
on 16 Feb 2018

As you can see from the Username, it's code I wrote by myself.

Maybe if you ask more clearly what you don't understand, I will be able to explain.

Best wishes

Torsten.

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

Start Hunting!