# Problem about mixture of ode-pde

15 views (last 30 days)
ZIYI LIU on 23 Jul 2020
Commented: Paul Maurerlehner on 24 Jan 2022
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
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
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.

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
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:
% 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
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.