matlab ode45求解齿轮动力学,Matlab拟合动力学参数遇到问题（ode45）

function [rho2,F] = rho2_F(k,r,ss,Ne,Np)

sy = sum(r.^2);

rho2 = 1 - ss./sy

F = (sy - ss).*(Ne-Np)./(Np.*ss) % F：F比

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')

fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))

fprintf('\tE1 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))

fprintf('\tk2 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))

fprintf('\tE2 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))

fprintf('\tk3 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))

fprintf('\tE3 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))

fprintf('\tk4 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))

fprintf('\tE4 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))

fprintf('\tk5 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9))

fprintf('\tE5 = %.4f ± %.4f\n',k(10),ci(10,2)-k(10))

fprintf('\tk6 = %.4f ± %.4f\n',k(11),ci(11,2)-k(11))

fprintf('\tE6 = %.4f ± %.4f\n',k(12),ci(12,2)-k(12))

fprintf('\tk7 = %.4f ± %.4f\n',k(13),ci(13,2)-k(13))

fprintf('\tE7 = %.4f ± %.4f\n',k(14),ci(14,2)-k(14))

fprintf('\tk8 = %.4f ± %.4f\n',k(15),ci(15,2)-k(15))

fprintf('\tE8 = %.4f ± %.4f\n',k(16),ci(16,2)-k(16))

fprintf('\tk9 = %.4f ± %.4f\n',k(17),ci(17,2)-k(17))

fprintf('\tE9 = %.4f ± %.4f\n',k(18),ci(18,2)-k(18))

function f = fy(k,x0,y0,ts,tt)

for i=1:42,%42组试验点

[t Xsim] = ode45(@g,[0,ts(i)],x0(i,:),[],k,tt(i));

i

ysim(i,1) = Xsim(end,1);

ysim(i,2) = Xsim(end,2);

ysim(i,3) = Xsim(end,3);

ysim(i,4) = Xsim(end,4);

ysim(i,5) = Xsim(end,5);

ysim(i,6) = Xsim(end,6);

ysim(i,7) = Xsim(end,7);

ysim(i,8) = Xsim(end,8);

ysim(i,9) = Xsim(end,9);

ysim(i,10) = Xsim(end,10);

ysim(i,11) = Xsim(end,11);

ysim(i,12) = Xsim(end,12);

%12个组分

end;

f = [ysim(:,:)-y0(:,:)];

function dCdt =g(t,C,k,tt)              % ODE模型

r1=k(1).*exp(-k(2)./8.314.*(1/tt-1/653)).*C(1).*C(2)-k(1).*exp(-k(2)./8.314.*(1/tt-1/653))./(0.001*tt*tt-1.514*tt+514.6).*C(4);

r2=k(3).*exp(-k(4)./8.314.*(1/tt-1/653)).*C(2).*C(2)-k(3).*exp(-k(4)./8.314.*(1/tt-1/653))./(0.00002*tt*tt-0.033*tt+11.44).*C(5);

r4=k(5).*exp(-k(6)./8.314.*(1/tt-1/653)).*C(2).*C(3)-k(5).*exp(-k(6)./8.314.*(1/tt-1/653))./(0.00003*tt*tt-0.046*tt+15.94).*C(6);

r5=k(7).*exp(-k(8)./8.314.*(1/tt-1/653)).*C(2).*C(4)-k(7).*exp(-k(8)./8.314.*(1/tt-1/653))./(0.00003*tt*tt-0.041*tt+14.14).*C(7);

r6=k(9).*exp(-k(10)./8.314.*(1/tt-1/653)).*C(3).*C(4)-k(9).*exp(-k(10)./8.314.*(1/tt-1/653))./(0.00003*tt*tt-0.034*tt+11.86).*C(8);

r8=k(11).*exp(-k(12)./8.314.*(1/tt-1/653)).*C(6);

r9=k(13).*exp(-k(14)./8.314.*(1/tt-1/653)).*C(7);

r10=k(15).*exp(-k(16)./8.314.*(1/tt-1/653)).*C(8) ;

r11=k(17).*exp(-k(18)./8.314.*(1/tt-1/653)).*C(2) ;

dC1dt=-r1;

dC2dt=-r1-2*r2-r4-r5-r11;

dC3dt=-r4-r6;

dC4dt=r1-r5-r6;

dC5dt=r2;

dC6dt=r4-r8;

dC7dt=r5-r9;

dC8dt=r6-r10;

dC9dt=r8;

dC10dt=r9;

dC11dt=r10;

dC12dt=r11;

dCdt = [dC1dt;dC2dt;dC3dt;dC4dt;dC5dt;dC6dt;dC7dt;dC8dt;dC9dt;dC10dt;dC11dt;dC12dt];