| 查看: 640 | 回复: 0 | ||
[求助]
求1stOpt程序代跑，求解微分方程组系数，可有偿，捉急捉急
|
|
需要拟合k1-k8的参数值 以下是我编写的matlabm文件: function QQQ_4_13_NJ_cs_2 format long clear all clc tspan = [0 0.25 1 3 9 22 35 91 135 ]; x0 = [19750 0 0 0 5250]; k0 = [0.147;0.19;0.03 ;0.02;0.01;0.000001184;3.7;13]; lb = [0;0;0;0;0;0;0;0]; ub = []; data=... [ 0.25 21300 100.225875 46.63575 16.168 3536.970375 1 22050 40.657785 14.123025 8.6688 2874.46475 3 21750 101.43456 48.06252 19.9383 3052.32539 9 19500 279.2715 190.106 23.3415 4908.7372 22 20200 390.3344 403.5 38.5308 3759.5435 35 18550 245.784 262.899 11.2602 5756.7424 91 17950 296.3025 363.12 20.31 6015.16875 135 17500 279.4176 467.418 11.8272 6891.531067 ]; yexp = data(:,2:6); [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk(1) = %.9f ± %.9f\n',k(1),ci(1,2)-k(1)) fprintf('\tk(2) = %.9f ± %.9f\n',k(2),ci(2,2)-k(2)) fprintf('\tk(3) = %.9f ± %.9f\n',k(3),ci(3,2)-k(3)) fprintf('\tk(4) = %.9f ± %.9f\n',k(4),ci(4,2)-k(4)) fprintf('\tk(5) = %.9f ± %.9f\n',k(5),ci(5,2)-k(5)) fprintf('\tk(6) = %.9f ± %.9f\n',k(6),ci(6,2)-k(6)) fprintf('\tk(7) = %.9f ± %.9f\n',k(7),ci(7,2)-k(7)) fprintf('\tk(8) = %.9f ± %.9f\n',k(8),ci(8,2)-k(8)) fprintf('The sum of the squares is:%.9e\n\n',resnorm) 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); Xsim4=Xsim(:,5); ysim(:,1) = Xsim1(2:end); ysim(:,2) = Xsim2(2:end); ysim(:,3) = Xsim3(2:end); ysim(:,4) = Xsim4(2:end); ysim(:,5) = Xsim4(2:end); size(ysim(:,1)); size(ysim(:,2)); size(ysim(:,3)); size(ysim(:,4)); size(ysim(:,5)); size(yexp(:,1)); size(yexp(:,2)); size(yexp(:,3)); size(yexp(:,4)); size(yexp(:,5)); f1=ysim(:,1)-yexp(:,1); f2=ysim(:,2)-yexp(:,2); f3=ysim(:,3)-yexp(:,3); f4=ysim(:,4)-yexp(:,4); f5=ysim(:,5)-yexp(:,5); f = [f1;f2;f3;f4;f5]; function dCdt = KineticsEqs(t,C,k) % ODE模型方程 dCAdt =-(k(3)+k(5)+k(7))*C(1)+k(1)*C(2)+k(2)*C(3)+k(2)*C(4)+k(8)*C(5); dCBdt = k(3)*C(1)-(k(1)+k(6))*C(2)+k(4)*C(3)+k(4)*C(4); dCCdt = k(5)*C(1)+k(6)*C(2)-(k(2)+k(4))*C(3); dCDdt = k(5)*C(1)+k(6)*C(2)-(k(2)+k(4))*C(4); dCEdt = k(7)*C(1)-k(8)*C(5); dCdt = [dCAdt;dCBdt;dCCdt;dCDdt;dCEdt]; |
回复此楼» 猜你喜欢
德国亥姆霍兹Hereon中心招收两位医用镁合金腐蚀与LPSO相变方向2026公派博士生
已经有0人回复
-
推荐一款可以AI辅助写作的Latex编辑器SmartLatexEditor，超级好用，推荐试试
已经有11人回复
-
物理学I论文润色/翻译怎么收费?
已经有276人回复
-
2026-CJ开始申报了
已经有1人回复
-
西安电子科学大学杭州研究院刘丽香教授招收智能多模态传感器和微型储能器件方向博士
已经有13人回复
-
重庆交通大学光子学微结构与器件课题组招收2026年硕士研究生信息
已经有1人回复
-
一志愿郑大材料学硕298分，求调剂
已经有6人回复
-
寻合作：应力腐蚀多尺度模拟
已经有3人回复
-
考研交流
已经有0人回复
5