MATLAB Answers

How can i solve a system coupled of PDE with an ODE using finite difference?

16 views (last 30 days)
Moji
Moji on 12 Feb 2018
Commented: Torsten on 16 Feb 2018
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)

Accepted Answer

Torsten
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.
  10 Comments
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!