MATLAB Answers

Problem about mixture of ode-pde

3 views (last 30 days)
ZIYI LIU
ZIYI LIU on 23 Jul 2020
Commented: ZIYI LIU on 8 Oct 2020
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

Show 1 older comment
ZIYI LIU
ZIYI LIU on 23 Jul 2020
Here is a system of two reactions, and 1, 2 means terminations (poles) of the system region. g and G means the same chemical, but G means inactive chemical and g means active chemical, and the same with c and C.
k+ and k-(kbar) are activation and deactivation rates of C(c), note that k+ depends on both the concentration of c at the poles and the amount of g at the poles. kon and koff are activation and deactivation rates of G(g), and kon depends on both the concentration of c at the poles. kcat descibes the relations between c and k+, and K describes the relations between c and kon. D means the diffusion rates of g and c.
Thanks!
Bill Greene
Bill Greene on 24 Jul 2020
Can you please provide values for all the parameters in the model?
ZIYI LIU
ZIYI LIU on 24 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

  7 Comments

Show 4 older comments
Bill Greene
Bill Greene on 29 Sep 2020
When you say, "I don't know how to modify the code", are you referring to the example script I show above? Use of pde1d assumes that you are reasonably familiar with pdepe. Specifically, what part of the pdepe or pde1d documentation are you confused about?
ZIYI LIU
ZIYI LIU on 6 Oct 2020
Hi Bill, I compared pdepe with pde1d (example script you show above), and I checked the pde1d code that you show above and I think they are right, such as the boundary conditions, and initial conditions...
However, I can't use pde1d get the right figures which I show above.
Yes, I am confused about some parts in your example code:
  1. c=sol(:,:,1); g=sol(:,:,2); Do these two terms describe C and G in the pde?
  2. plot(t, c(:,1)); What's the meaning of this term? Does it describe the concentration of C in the left end but not c1?
  3. How can I get figures of c1, c2, g1, g2 with time?
Thank you!
ZIYI LIU
ZIYI LIU on 8 Oct 2020
Hi Bill, Thank you so much! I got the figures! Thanks!
Best,
Ziyi Liu

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!