matlab求四元方程导数,matlab 求指点 动力学方程拟合过程中导数的获取 - 计算模拟 - 小木虫 - 学术 科研 互动社区...

本人matlab不通，勉强做了一个程序能够达到预期目的，可是老师想要中间运算的导数值，也就是dxdt的值，不知道如何填补表达式实现，望各位高手出手相助，不胜感激。

主程序：

format short

global Umax Ks Kp Y1 Y2 Uo Ko D Ux Kx Ku Ki;

Ki=0.0001;Ku=0.0006;

k0=[0.2 0.5 0.1 10 10 0.01 0.01 0.1 0.01 0.01]; %参数的初始值

x0=[0.0514 0.250 0]; %菌体浓度、底物浓度和产物浓度的初始值

t1=[0 6 12 18 24 30 36 42 48 54 60 66]';%发酵时间的实验数据

tspan=[0 6 12 18 24 30 36 42 48 54 60 66 ]';%时间间隔

%yexp实验数据[菌体浓度x1、底物浓度x2和产物浓度x3]

yexp=[[0.0514 0.0711 0.0997 0.1887 0.3478 0.6634 0.8711 1.1023 1.2336 1.4239 1.5849 1.6739];[0.25 0.2447 0.2366 0.2258 0.2196 0.203 0.1872 0.1698 0.1453 0.1124 0.0855 0.0608];[0 0.0012 0.0021 0.0039 0.0088 0.0137 0.0183 0.0268 0.0361 0.0438 0.0557 0.0643]]';

lb=[0 0 0 0 0 0.0216 0.001 0 0.001 0.001];ub=[1 10 10 50 50 2 1 1 0.1 0.1];%参数的下、上限

[k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@ObjFunc4LNL103,k0,lb,ub,[],x0,yexp);%非线性拟合

y1=[yexp(:,1)];y2=[yexp(:,2)];y3=[yexp(:,3)]';

[t4plot,x4plot]=ode45(@kineticseqs103,[0 100],x0,[],k);

plot(t1,y1,'b*',t4plot,x4plot(:,1),'k-'),xlabel('T(h)'),ylabel('Cx/(g/L)');

figure

plot(t1,y2,'g*',t4plot,x4plot(:,2),'k-'),xlabel('T(h)'),ylabel('Cs/(mmol/L)');

figure

plot(t1,y3,'r*',t4plot,x4plot(:,3),'k-'),xlabel('T(h)'),ylabel('Cp/(mmol/L)');

disp(k)

子程序：

1)

function f=ObjFunc4LNL103(k,x0,yexp)

tspan=[0 6 12 18 24 30 36 42 48 54 60 66];

[t1,x]=ode45(@kineticseqs103,tspan,x0,[],k);

y(:,1)=x(:,1);y(:,2)=x(:,2);y(:,3)=x(:,3);

f1=y(: ,1)-yexp(: ,1);f2=y(: ,2)-yexp(: ,2);f3=y(: ,3)-yexp(: ,3);

f=[f1*1 f2*5 f3*10];

2)

function dxdt=kineticseqs103(t,x,k) %模型方程

global Umax Ks Kp Y1 Y2 Uo Ko D Ux Kx Ku Ki

Umax=k(1); Ks=k(2); Kp=k(3); Y1=k(4); Y2=k(5); Ux=k(6); Kx=k(7); D=k(8); Uo=k(9); Ko=k(10);

dxdt1=Umax*x(2)*x(1)*(1-Kp*x(3))/(Ks+x(2));

if dxdt1<0

dxdt1=0;

end

dxdt2=-1/Y1*dxdt1-Ux*x(2)*x(1)/(Kx+x(2))*(Ki/(Ki+x(3)));

dxdt3=D*(-dxdt2-1/Y2*dxdt1-Uo*x(2)*x(1)/(Ko+x(2))*(Ku/(Ku+x(3))));

if dxdt3<0

dxdt3=0;

end

dxdt=[dxdt1;dxdt2;dxdt3];

只要能显示对应时间t1,这些点的dxdt值就可以了。谢谢