| 查看: 765 | 回复: 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]; |
回复此楼» 猜你喜欢
求助中医院临床中医中标的青年标书学习，谢谢
已经有6人回复
-
浙大物院表面物理课题组长期招收普博及直博学生
已经有259人回复
-
物理学I论文润色/翻译怎么收费?
已经有146人回复
-
【征稿】Springer Nature热电专题 Discover Materials / Discover Chemistry
已经有47人回复
-
25届硕士申请，材料计算方向
已经有2人回复
-
课题组招1-2名硕士
已经有1人回复
-
机时越来越贵之后，我把一部分 DFT 流程搬到了 AI 计算平台上
已经有0人回复
-
斯坦福大学与Elsevier世界前2%学者新排名查询
已经有1人回复
-
27申博推荐：南京航空航天大学国际前沿院光学方向招收博士生！
已经有19人回复
-
华南师范大学+双一流学科物理学+广东省杰青团队2027年博士招生+机器学习驱动材料设计
已经有1人回复
-
apl投稿状态
已经有2人回复
20