Problem about mixture of ode-pde

8 次查看（过去 30 天）
ZIYI LIU
ZIYI LIU 2020-7-23
评论： Paul Maurerlehner 2022-1-24
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 个评论
显示 2更早的评论
Bill Greene
Bill Greene 2020-7-24
Can you please provide values for all the parameters in the model?
ZIYI LIU
ZIYI LIU 2020-7-24
编辑：ZIYI LIU 2020-7-25
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bill Greene
Bill Greene 2020-7-29
编辑：Bill Greene 2020-10-8
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 个评论
显示 9更早的评论
Bill Greene
Bill Greene 2022-1-24
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
Paul Maurerlehner 2022-1-24
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.

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by