9 views (last 30 days)

Show older comments

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
on 29 Jul 2020

Edited: Bill Greene
on 8 Oct 2020

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

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

Start Hunting!