| 查看: 686 | 回复: 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 |
回复此楼» 猜你喜欢
北伊利诺伊大学(Northern Illinois University)机械系招收全奖博士生
已经有0人回复
-
第五届粉末衍射晶体结构解析及Rietveld精修高级研讨会
已经有45人回复
-
物理学I论文润色/翻译怎么收费?
已经有132人回复
-
FDTD软件里的时域图像画图
已经有2人回复
-
复杂量子系统的新物态和新机理讲习班
已经有0人回复
-
感谢小木虫，已完成目标
已经有0人回复
-
零点能
已经有2人回复
-
科研求助
已经有0人回复
-
【25申博】博导推荐求助
已经有6人回复
-
青岛大学李山东教授招生
已经有2人回复
独孤神宇
版主 (知名作家)
- 应助: 489 (硕士)
- 贵宾: 0.008
- 金币: 29623.8
- 散金: 802
- 红花: 121
- 沙发: 1
- 帖子: 5358
- 在线: 850.4小时
- 虫号: 3522474
- 注册: 2014-11-06
- 性别: GG
- 专业: 机械动力学
- 管辖: 计算模拟
2楼2021-01-05 21:03:26
3楼2021-01-05 21:07:08
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15356.8
- 散金: 101
- 红花: 234
- 帖子: 3409
- 在线: 1222.7小时
- 虫号: 291104
- 注册: 2006-10-28
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
自制奶酪: 金币+20, ★有帮助 2021-01-05 21:39:36
独孤神宇: 金币+5, 鼓励交流 2021-01-06 17:08:32
感谢参与，应助指数 +1
自制奶酪: 金币+20, ★有帮助 2021-01-05 21:39:36
独孤神宇: 金币+5, 鼓励交流 2021-01-06 17:08:32
|
参考下1stOpt的计算结果： Weighted Root of Mean Square Error (RMSE): 3.52215359203778 Weighted Sum of Squared Residual: 434.194807406662 Correlation Coef. (R): 0.967264353331916 R-Square: 0.935600329226609 Determination Coef. (DC): 0.698480780131114 F-Statistic: -0.866225958773538 Parameter Best Estimate -------------------- ------------- k1 1.18206935637483 k2 1.23804230026633 k3 0.392064236557302 k4 0.237399839174299 k5 2.4090372987211E-18 k6 1.44475738382174E-24 k7 0.719127412563461 f Initial Value 19.2237430015159 |
4楼2021-01-05 21:34:03
5楼2022-11-17 13:58:31
5