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

7 次查看（过去 30 天）
Moji
Moji 2018-2-12
评论： Kuldeep Malik 2023-8-10
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
Torsten 2018-2-13
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 个评论
显示 8更早的评论
Torsten
Torsten 2018-2-16
Your new code is too much different from the one I linked to.
Sorry, but I don't have the time to dive in that deep.
Best wishes
Torsten.
Torsten
Torsten 2018-2-16
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by