![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/336220/image.png)
Problem about mixture of ode-pde
조회 수: 8 (최근 30일)
이전 댓글 표시
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
채택된 답변
Bill Greene
2020년 7월 29일
편집: Bill Greene
2020년 10월 8일
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
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
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에서 Boundary Conditions에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!