常微分方程组求解
2017-05-25
2017
计算得到的了反应速率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，非R2Comments