| 查看: 667 | 回复: 4 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
用matlab求动力学参数已有2人参与
|
||
|
我不知道MATLAB怎么输入求动力学函数的方程，抄了木虫大佬的内容，请大家帮我把它弄到能运行出来，动力学参数拟合问题，求助大家~~~ function odes_fit format long clear all clc k0=[0 0 0 0 0 0 0 0];%参数初值 lb=[0 0 0 0 0 0 0 0];ub=[+inf +inf +inf +inf +inf +inf +inf +inf];%lb、ub:参数下限和上限 x0=[0 0 0 0 0]; data=... [0.5 1 1.5 2 3 4 5 6 18.7 6.5 2.7 2.1 1.7 1.5 1.2 0.2; 48.9 31.3 15.3 6.9 2.2 1.8 1.7 0.5; 19.3 34.7 41.2 40.9 37.2 31.7 26.1 19.4; 9.7 19.6 28.7 34.5 38.4 38.2 36.3 33.4; 0.2 2.6 5.7 9.3 17.1 24.8 30.5 36.7; ]; x0=data(1,2:end); tspan = [data(:,1)']; yexp = [data(2:end,2) data(2:end,3) data(2:end,4) data(2:end,5) data(2:end,6) data(2:end,7) data(2:end,8)]; %使用函数fmincon()进行参数估计 [k,fval,flag]=fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.9f \n',k(1)) fprintf('\tk2 = %.9f \n',k(2)) fprintf('\tk3 = %.9f \n',k(3)) fprintf('\tk4 = %.9f \n',k(4)) fprintf('\tk5 = %.9f \n',k(5)) fprintf('\tk6 = %.9f \n',k(6)) fprintf('The sum of the squares is:%.le\n\n',fval) k_fmincon=k; %使用函数lsqnonlin()进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian]=... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp); ci=nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值:\n'),output %以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计 k0=k_fmincon; [k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp); ci=nlparci(k,residual,jacobian); fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值:\n') output figure(1) ts=0  (max(tspan)-min(tspan))/100):max(tspan);[ts, ys] = ode45(@KineticsEqs,ts,x0,[],k); yy = [data(:,2) data(:,3) data(:,4) data(:,5) data(:,6)]; figure(1) plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo'); figure(2) plot(ts,ys(:,2),'r',tspan,yy(:,2),'ro'); figure(3) plot(ts,ys(:,3),'k',tspan,yy(:,3),'ko'); figure(4) plot(ts,ys(:,4),'g',tspan,yy(:,4),'ko'); figure(5) plot(ts,ys(:,4),'g',tspan,yy(:,4),'ko'); figure(6) plot(ts,ys(:,6),'m',tspan,yy(:,6),'mo'); %legend('C1的计算值','C1的实验值','C2的计算值','C2的实验值','C3的计算值','C3的实验值','C4的计算值','C4的实验值','C5的计算值','C5的实验值','C6的计算值','Location','best'); function f = ObjFunc(k,tspan,x0,yexp) % 目标函数 [t, Xsim] = ode45(@KineticsEqs,tspan,x0,[],k); Xsim1=Xsim(:,1); Xsim2=Xsim(:,2); Xsim3=Xsim(:,3); Xsim4=Xsim(:,4); Xsim5=Xsim(:,5); Xsim6=Xsim(:,6); ysim(:,1) = Xsim1(2:end); ysim(:,2) = Xsim2(2:end); ysim(:,3) = Xsim3(2:end); ysim(:,4) = Xsim4(2:end); ysim(:,5) = Xsim5(2:end); ysim(:,6) = Xsim6(2:end); f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3)) (ysim(:,5)-yexp(:,5)) ... (ysim(:,6)-yexp(:,6))]; function dCdt = KineticsEqs(t,C,k) % ODE模型方程 C1=C(1);C2=C(2);C3=C(3);C4=C(4);C5=C(5);C6=C(6); k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);k7=k(7); dC1dt = -k1*C1-k7*C1; dC2dt = k1*C1-k2*C2-k5*C2; dC3dt = k2*C2-k3*C3-k6*C3; dC4dt = k3*C3-k4*C5; dC5dt = k4*C5; dC6dt = k5*C2+k6*C3+k7*C1; dCdt = [dC1dt;dC2dt;dC3dt;dC4dt;dC5dt;dC6dt]; @beefly |
回复此楼» 猜你喜欢
法国博士招聘
已经有0人回复
-
法国博士招聘
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有145人回复
-
法国博士招聘
已经有0人回复
-
法国博士招聘
已经有0人回复
-
超导三元凸包图
已经有0人回复
-
最近遇到这样一个问题
已经有10人回复
-
论文接收，散金币了！
已经有176人回复
-
vasp算居里温度
已经有0人回复
-
郑州大学光电器件方向招收一名博士
已经有6人回复
3楼2021-01-05 21:07:08
100