| 查看: 1482 | 回复: 2 | ||
[求助]
Matlab拟合反应动力学参数结果偏差很大啊，求大神指点程序当如何修改 已有1人参与
|
|
最近一直在拟合反应动力学参数，花了很大功夫试着编好了代码，却始终不能得出文献中的参数结果，求大神帮忙指点一下程序。在此拜谢了。 code: function KineticsEst5 clear all clc k0 = [0.5 0.5 0.5 0.5 0.5 0.5]; % 参数初值 lb = [0 0 0 0 0 0]; % 参数下限 ub = [+inf +inf +inf +inf +inf +inf]; % 参数上限 x0 = [4.96 24.43 26.32 17 38.72]; ExpData = ... [ 0 4.96 24.43 26.32 5 5.635 28.29 30.715 10 6.31 32.15 35.11 15 6.77 34.225 34.98 20 7.23 36.3 34.85 25 7.065 39.285 33.525 30 6.9 42.27 32.2 35 7.08 45.46 29.7 40 7.26 48.65 27.2 50 8 51.215 24.495 60 8.74 53.78 21.79 75 8.605 58.01 19.12 90 8.47 62.24 16.45 ] yexp = ExpData(:,2:4); % 使用函数fmincon()进行参数估计 [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp); fprintf(\\\'\\\\n使用函数fmincon()估计得到的参数值为:\\\\n\\\') fprintf(\\\'\\\\tk1 = %.4f\\\\n\\\',k(1)) fprintf(\\\'\\\\tk2 = %.4f\\\\n\\\',k(2)) fprintf(\\\'\\\\tk3 = %.4f\\\\n\\\',k(3)) fprintf(\\\'\\\\tk4 = %.4f\\\\n\\\',k(4)) fprintf(\\\'\\\\tk5 = %.4f\\\\n\\\',k(5)) fprintf(\\\'\\\\tk6 = %.4f\\\\n\\\',k(6)) fprintf(\\\' The sum of the squares is: %.1e\\\\n\\\\n\\\',fval) k_fmincon = k; % 以函数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\\\') fprintf(\\\'\\\\tk1 = %.4f ± %.4f\\\\n\\\',k(1),ci(1,2)-k(1)) fprintf(\\\'\\\\tk2 = %.4f ± %.4f\\\\n\\\',k(2),ci(2,2)-k(2)) fprintf(\\\'\\\\tk3 = %.4f ± %.4f\\\\n\\\',k(3),ci(3,2)-k(3)) fprintf(\\\'\\\\tk4 = %.4f ± %.4f\\\\n\\\',k(4),ci(4,2)-k(4)) fprintf(\\\'\\\\tk5 = %.4f ± %.4f\\\\n\\\',k(5),ci(5,2)-k(5)) fprintf(\\\'\\\\tk6 = %.4f ± %.4f\\\\n\\\',k(6),ci(6,2)-k(6)) fprintf(\\\' The sum of the squares is: %.1e\\\\n\\\\n\\\',resnorm) % ------------------------------------------------------------------ function f = ObjFunc4Fmincon(k,x0,yexp) tspan = [0,5,10,15,20,25,30,35,40,50,60,75,90]; [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1:3) = x(:,1:3); f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2); % ------------------------------------------------------------------ function f = ObjFunc4LNL(k,x0,yexp) tspan = [0,5,10,15,20,25,30,35,40,50,60,75,90]; [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1:3) = x(:,1:3); f1 = y(:,1) - yexp(:,1); f2 = y(:,2) - yexp(:,2); f3 = y(:,3) - yexp(:,3); f = [f1; f2; f3]; % ------------------------------------------------------------------ function dxdt = KineticEqs(t,x,k) dxdt = ... [ (k(1)*x(1)+k(6)*x(3)) (k(2)*x(1)+k(5)*x(3)) (k(3)*x(1)+k(4)*x(2)-(k(5)+k(6))*x(3)) (-(k(1)+k(2)+k(3))*x(1)) (-k(4)*x(2)) ]; |
回复此楼» 猜你喜欢
散金币祈福
已经有186人回复
-
求助，如何提取ELK的rt-TDDFT在某一时刻的自旋密度分布和电子密度分布
已经有1人回复
-
物理学I论文润色/翻译怎么收费?
已经有59人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
深圳大学应用技术学院招聘凝聚态物理博士后
已经有112人回复
-
深圳大学2027级光电信息工程/物理学研究生推免
已经有69人回复
-
求助中医院临床中医中标的青年标书学习，谢谢
已经有4人回复
-
浙大物院表面物理课题组长期招收普博及直博学生
已经有257人回复
-
【征稿】Springer Nature热电专题 Discover Materials / Discover Chemistry
已经有39人回复
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
2楼2015-10-01 19:12:13
3楼2015-10-01 20:03:33
