把梦想照进现实

常微分方程组求解

2017.05.25

计算得到的了反应速率k,需要得到常微分方程组的解。在网上查找了下,整理如下:
新建rigid.m函数

function dy=rigid(t,y)
k1=0.390161;
k2=0.115611;
k3=0.016353;
k4=0.0562543;
k5=0.281573;
k6=0.4055;
k7=0.0172644;
k8=5.86848e-6;
k9=1.31968e-9;
n=0.407151;
dy=zeros(4,1);
dy(1)=k7*y(2)*n+k9*y(4)+k8*y(3);
dy(2)=-k1*y(2)*n-k3*y(2)*n-k7*y(2)*n+k2*y(3)+k4*y(3);
dy(3)=-k4*y(3)-k8*y(3)-k6*y(3)+k3*y(2)*n+k5*y(4);
dy(4)=-k2*y(4)-k9*y(4)-k5*y(4)+k1*y(2)*n+k6*y(3);

保存后,运行。

[T,Y]=ode45(‘rigid’,[0 90],[0 100 0 0]);
plot(T,Y(:,1),’-‘,T,Y(:,2),’*’,T,Y(:,3),’+’,T,Y(:,4),’.’)

为了得到拟合曲线和原始数字的方差

%差值取数
tt=[3;15;30;60;90];               %试验记录时间点
s3=interp1(T,Y,tt)                %该试验时间点下的拟合值,自动导入
s33=[8.097801092 8.52400115 6.748167577 6.748167577 9.234334579];                 %原始数字
c = corrcoef(s3(:,1),s33)         %相关系数R,非R2
Comments
Write a Comment