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

16 views (last 30 days)
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)

Torsten on 13 Feb 2018
Take a look at the answer provided here:
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.