| 查看: 658 | 回复: 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]; |
回复此楼» 猜你喜欢
欢迎加入课题组
已经有0人回复
-
散金币，求好运，祝面上顺利！
已经有80人回复
-
物理学I论文润色/翻译怎么收费?
已经有171人回复
-
华东师范大学芯片设计徐珑真实水平咋样
已经有1人回复
-
调剂
已经有0人回复
-
【新加坡】纳米电子器件项目组有“联合培养博士生”名额
已经有0人回复
-
中德博士后2026项目5月30日截止- 国内博士后入站、德国两年全职
已经有0人回复
-
面向全球申请人的中德博士后2026项目5月30日截止- 国内博后入站德国两年全职
已经有0人回复
-
【新加坡】纳米电子器件项目组有“联合培养博士生”名额
已经有2人回复
-
重庆交大材料学院26年硕士招生复试实施细则已出。欢迎加入光子学微结构与器件课题组。
已经有0人回复
-
北工大宋晓艳教授团队招收2026年9月入学博士【学博和专博名额均有】！
已经有1人回复
40