Problem about mixture of ode-pde

16 views (last 30 days)
Hello,
I have a system that I want to solve numerically (attached is the pde-ode file). I know a little about how to solve ode and pde separately, but I don't know how to combine ode and pde parts in MATLAB.
Here is my code: (I am sorry for my rough code)
T=120;
c10 = 1.3;
c20 = 0.03;
g10 = 0.07;
g20 = 1.37;
[t, state_variable]=ode45(@LV,[0 T],[c10,c20,g10,g20]);
c1 = state_variable(:,1);
c2 = state_variable(:,2);
g1 = state_variable(:,3);
g2 = state_variable(:,4);
%[c1,t]
x = linspace(0,L,20);
%t = [linspace(0,0.05,20), linspace(0.5,5,10)];
m = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);
hold on
plot(t,c1,'LineWidth',2)
plot(t,c2,'LineWidth',2)
% Various commands to make graph clearer
h=legend('c1','c2','g1','g2');
xlabel('Time','Fontsize',16)
ylabel('concentration','Fontsize',16)
set(gca,'XTick')
set(h,'Fontsize',20)
function f=LV(t,state_variable)
c1=state_variable(1);
c2=state_variable(2);
g1=state_variable(3);
g2=state_variable(4);
%here i deleted many parameters
% Equations
dc1 =(k0+kcat*c1^n)*g1*C(0,t)-kbar*c1;
dc2 =(k0+kcat*c2^n)*g2*C(L,t)-kbar*c2;
dg1=kon/(1+K*c1^m)*G(0,t)-koff*g1;
dg2=kon/(1+K*c2^m)*G(L,t)-koff*g2;
%f=[dx;dy;];
f=[dc1;dc2;dg1;dg2;];
end
function [c,f,s] = heatpde(x,t,C,dCdx)
c = 1/Dc;
f = dCdx;
s = 0;
end
function C0 = heatic(x)
C0 = 0.17;
end
function [pl,ql,pr,qr] = heatbc(xl,Cl,xr,Cr,t)
pl = -Dc*Cl+dc1;
ql = 0;
pr = -Dc*Cr-dc2;
qr = 0;
end
Thank you very much!
Best,
Ziyi
  4 Comments
ZIYI LIU
ZIYI LIU on 24 Jul 2020
Edited: ZIYI LIU on 25 Jul 2020
Below are my values in the code I wrote by using MOL just now。 But I cannot let C,C,g,G change with time and position, my result is just the initial result. Maybe I didn't make correct connections between every part. Thanks!
%main program
clear all
clc
close all
% Data
L=1;
Nz=100; %n=Nz+1
Dc=3;
Dg=3;
dz=L/Nz;
Ctot=1.58;
Gtot=1.5;
c10=1.3;
c20=0.03;
g10=0.07;
g20=1.37;
k0 =0.1;
kcat= 40;
kbar = 1;
m =2;
n=2;
kon=1;
K=8;
koff=0.9;
%g1=1;%
%g2=0.1;%
t = [0:0.01:50];
z= [0:0.01:L];%?I am not sure about linspace of t and (0,L,Nz)?
%initial condition at t=0
IC=zeros(Nz+1,1);%or zeros(Nz+1?1)?zeros(1,Nz)? and should if be (Nz+1,1)like the vedio
%IC(1:Nz+1)=(Ctot-c10-c20)/Nz;% i assume.or(1:Nz+1)
IC(1:Nz+1)=0.25;% i assume.or(1:Nz+1)
IG=zeros(Nz+1,1);
%IG(1:Nz+1)=(Gtot-g10-g20)/Nz;
IG(1:Nz+1)=0.06;% i assume.or(1:Nz+1)
Ic=zeros(Nz+1,1);%or(Nz+1,1).
Ic(1)=c10;
Ic(Nz+1)=c20;
Ig=zeros(Nz+1,1);%or(Nz+1,1).
Ig(1)=g10;
Ig(Nz+1)=g20;
Iy=[IC;Ic;IG;Ig];
%Solver
%[t, state_variable]=ode45(@LV,tspan,[c10,c20]);
%c1 = state_variable(:,1);
%c2 = state_variable(:,2);
[T ,Y]=ode15s(@(t,y) fun(t,y,z,Nz,Dc,Dg,dz,c10,c20,g10,g20,L,k0,kcat,kbar,kon,K,koff,n,m),t,Iy);%IC,[],Nz,Dc,dc1dt,dc2dt,c1,c2);
%plot 2D
hold on
plot(T, Y)%c1
hold on
%plot(T, Y(:,2*Nz+2))%c2
%plot(T, Y)%why only two lines?
ylabel('c1');
xlabel('time span');
%bar image 3D
%imagesc(z,t, Y(1:Nz+1))
%imagesc(z,t, Y(Nz+2:2*Nz+2))% z-xaxis,t-yaxis, and i only want c
%imagesc(t, Y(Nz+2:2*Nz+2)), i cannot get result(error)
%grid on
%xlabel('z position');
%ylabel('time span');
%colormap jet
%colorbar
%dispersion
function DyDt=fun(t,y,z,Nz,Dc,Dg,dz,c10,c20,g10,g20,L,k0,kcat,kbar,kon,K,koff,n,m)
C=zeros(Nz+1,1);
c=zeros(Nz+1,1);
G=zeros(Nz+1,1);
g=zeros(Nz+1,1);
%C=zeros(1,Nz+1);
%c=zeros(1,Nz+1);
DCDt=zeros(Nz+1,1);
DcDt=zeros(Nz+1,1);
DGDt=zeros(Nz+1,1);
DgDt=zeros(Nz+1,1);
DyDt=zeros(4*(Nz+1),1);
%zhalf=zeros(Nz,1);% WHY NEED THIS?
%D2CDz2=zeros(Nz-1,1);
C=y(1:Nz+1);
c=y(Nz+2:2*(Nz+1));
G=y(2*(Nz+1)+1:3*(Nz+1));
g=y(3*(Nz+1)+1:4*(Nz+1));
%zhalf(1:Nz)=(z(1:Nz)+z(2:Nz+1))/2;
%boundary conditions
C(1)=1/3*(4*C(2)-C(3)-2*dz/Dc*DcDt(1));
C(Nz+1)=1/3*(-2*dz/Dc*DcDt(Nz+1)+4*C(Nz)-C(Nz-1));
G(1)=1/3*(4*G(2)-G(3)-2*dz/Dg*DgDt(1));
G(Nz+1)=1/3*(-2*dz/Dg*DgDt(Nz+1)+4*G(Nz)-G(Nz-1));
for i=2:Nz
D2CDz2(i)=1/(dz^2)*(C(i+1)-2*C(i)+C(i-1));
D2GDz2(i)=1/(dz^2)*(G(i+1)-2*G(i)+G(i-1));
DcDt(i)=0;%as there is no c between 2:Nz, c is only existing in two poles
DgDt(i)=0;
end
%how to get i=Nz+1?=1?:%I am really not sure about this: (second order derivative)
D2CDz2(Nz+1)=1/(dz^2)*(-C(Nz-1)+2*C(Nz)-C(Nz+1));
D2CDz2(1)=1/(dz^2)*(C(3)-2*C(2)+C(1));
D2GDz2(Nz+1)=1/(dz^2)*(-G(Nz-1)+2*G(Nz)-G(Nz+1));
D2GDz2(1)=1/(dz^2)*(G(3)-2*G(2)+G(1));
%for time, (1:Nz+1)
%the rate of translation between C and c(G and g) at two poles.
DcDt(1) =(k0+kcat*c(1)^n)*g(1)*C(1)-kbar*c(1);
DcDt(Nz+1) =(k0+kcat*c(Nz+1)^n)*g(Nz+1)*C(Nz+1)-kbar*c(Nz+1);
DgDt(1)=kon/(1+K*c(1)^m)*G(1)-koff*g(1);
DgDt(Nz+1)=kon/(1+K*c(Nz+1)^m)*G(Nz+1)-koff*g(Nz+1);
for i=1:Nz+1
DCDt(i)=Dc*D2CDz2(i);
DGDt(i)=Dg*D2GDz2(i);
end
DyDt=[DCDt;DcDt;DGDt;DgDt];
%plot(t,c(1))
end
%problem is:
%I got from the workspace that my result is almost the initial conditions.
%C G c g aren't change.

Sign in to comment.

Accepted Answer

Bill Greene
Bill Greene on 29 Jul 2020
Edited: Bill Greene on 8 Oct 2020
I have developed a PDE solver, pde1dM, that Ibelieve can solve your coupled PDE/ODE system.
The solver runs in MATLAB and is similar to the standard pdepe solver but it allows a set of ODE to
be coupled to the set of PDE.
I have used your pdf document and example code to solve a problem which I think
is close to what you want to solve. I have attached this MATLAB script below.
Unfortunately, I don't understand your problem well enough to know if I have
accurately reproduced your intentions.
If you want to try or modify this example yourself, you can easily download pde1dM at this location and
run the example code shown below:
function matlabAnswers_7_25_2020
% Data
L=1;
Nz=100; %n=Nz+1
Dc=3;
Dg=3;
dz=L/Nz;
Ctot=1.58;
Gtot=1.5;
c10=1.3;
c20=0.03;
g10=0.07;
g20=1.37;
k0 =0.1;
kcat= 40;
kbar = 1;
m =2;
n=2;
kon=1;
K=8;
koff=0.9;
%g1=1;%
%g2=0.1;%
tFinal=15;
t=linspace(0,tFinal,30);
z=linspace(0,L,20);
xOde = [0 L]'; % couple ODE at the two ends
mGeom = 0;
%% pde1dM solver
pdef = @(z,t,u,DuDx) pdefun(z,t,u,DuDx,Dc,Dg);
ic = @(x) icfun(x);
odeF = @(t,v,vdot,x,u,DuDx,f, dudt, du2dxdt) ...
odeFunc(t,v,vdot,x,u,DuDx,f, dudt, du2dxdt, ...
k0, kcat, n, kbar, kon, K, m, koff);
odeIcF = @() odeIcFunc(c10,c20,g10,g20);
%figure; plot(x, ic(x)); grid; return;
[sol, odeSol] = pde1dM(mGeom,pdef,ic,@bcfun,z,t,odeF, odeIcF,xOde);
C=sol(:,:,1);
G=sol(:,:,2);
figure; plot(z, C(end,:)); grid;
title 'C at final time';
figure; plot(z, G(end,:)); grid;
title 'G at final time';
figure; plot(t, C(:,1)); grid;
title 'C at left end as a function of time';
figure; plot(t, C(:,end)); grid;
title 'C at right end as a function of time';
figure; plot(t, G(:,1)); grid;
title 'G at left end as a function of time';
figure; plot(t, G(:,end)); grid;
title 'G at right end as a function of time';
% plot ode variables as a function of time
figure;
hold on;
for i=1:4
plot(t, odeSol(:,i));
end
legend('c1', 'c2', 'g1', 'g2');
grid; xlabel('Time');
title 'ODE Variables As Functions of Time'
end
function [c,f,s] = pdefun(z,t,u,DuDx,Dc,Dg)
c = [1 1];
f = [Dc Dg]'.*DuDx;
s = [0 0]';
end
function c0 = icfun(x)
c0 = [.25 .06]';
end
function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t,v,vdot)
dc1dt = vdot(1);
dc2dt = vdot(2);
dg1dt = vdot(3);
dg2dt = vdot(4);
pl = [-dc1dt -dg1dt]';
ql = [1 1]';
pr = [dc2dt dg2dt]';
qr = [1 1]';
end
function R=odeFunc(t,v,vdot,x,u,DuDx,f, dudt, du2dxdt, ...
k0, kcat, n, kbar, kon, K, m, koff)
C1 = u(1,1);
C2 = u(1,2);
G1 = u(2,1);
G2 = u(2,2);
c1 = v(1);
c2 = v(2);
g1 = v(3);
g2 = v(4);
Dc1Dt =(k0+kcat*c1^n)*g1*C1-kbar*c1;
Dc2Dt =(k0+kcat*c2^n)*g2*C2-kbar*c2;
Dg1Dt=kon/(1+K*c1^m)*G1-koff*g1;
Dg2Dt=kon/(1+K*c2^m)*G2-koff*g2;
R=vdot-[Dc1Dt, Dc2Dt, Dg1Dt, Dg2Dt]';
end
function v0=odeIcFunc(c10,c20,g10,g20)
v0=[c10,c20,g10,g20]';
end
  11 Comments
Paul Maurerlehner
Paul Maurerlehner on 24 Jan 2022
You are right, at least the second equation at equations (5) is a PDE since it contains a spatial deriavative (of u_1). I could define all equations as PDEs, but I have the problem that I don't have boundary conditions for the auxiliary variables..
I postet it here, because this specific issue concerns the pde1dM solver which is not a standard Matlab function and the probleme here is very similar to mine.

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!