Problem about mixture of ode-pde

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

[Bill Greene]

I do not believe it is possible to solve this system with a combination of pdepe and ode45.

However, if you include the details you have omitted, e.g. definitions for

Dg, Dc, kbar, kon, koff, khat, ...

I would be interested in looking into it further. It is also not clear, in some places, how

your code coresponds to your mathematical description, e.g. .

[ZIYI LIU]

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]

Can you please provide values for all the parameters in the model?

[ZIYI LIU]

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.

Answer 1

[Bill Greene]

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

[ZIYI LIU]

Hi Bill, thank you so much for your answer, and I will try this soon!

[ZIYI LIU]

Hi Bill, I found some problems in the code. The initial position in the pictures doesn't consist with the initial number (c10=0.13,c20=003), and I don't understand the term,

function c0 = icfun(x)

c0 = [.1 1]';

end

The problem is that number in c0 = [.1 1] is unactive C (described as C) concentration in the cell, but not active C (described as c), and c0 = [.1 1] decides the initial position in the figures. I want to know how you add initial conditions in the code, such as c10 c20 g10 g20, and Ctot (C total) and Gtot are constant in the system.

The figures I want is that c1 and c2 are anti-correlated and oscillatory. Hope you can tell me how to deal with this problem.Thank you!

[Bill Greene]

The purpose of the icfun function is to set the initial values of the PDE dependent

variables-- in your equations, C and D. This function is defined exactly like the

corresponding pdepe function. As far as the actual values in my sample code, I

don't recall how I selected them.

You need to understand that my sample code is merely an attempt to get you started

in defining the appropriate code for your specific problem. You will still need to study

and understand the documentation for both pdepe and pde1d and make changes

to the sample code as you deem necessary.

Regarding the pde1d documentation, I am interested in improving that so if you find

places where it is unclear, please let me know and I will try to improve them.

[ZIYI LIU]

Hi Bill,

I checked the code and I think its good for my model, however, the results are different from what I want, and I don't know how to modify the code.

I am not sure about the meanings of c and g in the plot part. Are they small c and small g in my model? But the initial value of the plots doesn't consist with the initial value of c and g. Also the Ctot and Gtot are constant.

The result I want is in the picture below.

Also, I want to use the FEM pde toolbox but the boundary conditions are odes, and its difficult for me to use this. Can I use FEM (pde toolbox) to solve my model?

Thank you very much!

Best,

Ziyi Liu

[Bill Greene]

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]

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]

Hi Bill, Thank you so much! I got the figures! Thanks!

Best,

Ziyi Liu

[ZIYI LIU]

Hi Bill, I still have a question about the code.

I want to add some random distribution to the ode, but there is no result if I add rand(1).

Dc2Dt =1*(0.3*k0+1*kcat*c2^n)*g2*C2-1*kbar*c2+rand(1);

Could you help me with this?

Thank you very much!

Ziyi Liu

[Paul Maurerlehner]

Hi Bill,

Is there a way to couple the ODEs in the PDE function not only via the boundary conditions with your developed solver?

Because if I understand it right, you don't have the variables v and vdot available in the PDE function pdefun and in my problem (also see atteched pdf) I would require them.

Kind regards,

Paul

[Bill Greene]

The ODE variables are available to the PDE definition function. Look at section 3.2.3 in this pde1dM user manual that shows the function signature for the PDE function when there are ODE variables.

But looking at equations (5) in your doc, they look like PDE to me because they contain the u_2 and u_1 which are functions of x.

Why don't you post a new question here concerning your problem?

[Paul Maurerlehner]

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.