24小时热门版块排行榜    

查看: 1217  |  回复: 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
回复此楼
简单快乐
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

自制奶酪

铜虫 (小有名气)

引用回帖:
2楼: Originally posted by 独孤神宇 at 2021-01-05 21:03:26
这个代码不完整

lsqnonlin 函数 ObjFunc4LNL 缺失

您好,这个要怎么修改啊,求助
简单快乐
3楼2021-01-05 21:07:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 5 个回答

独孤神宇

版主 (知名作家)

【答案】应助回帖

感谢参与,应助指数 +1
这个代码不完整

lsqnonlin 函数 ObjFunc4LNL 缺失
数值计算
2楼2021-01-05 21:03:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与,应助指数 +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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

流浪在地球

新虫 (初入文坛)

5楼2022-11-17 13:58:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 东华理工大学化材专业26届硕士博士申请 +6 zlingli 2026-03-13 6/300 2026-03-15 20:00 by ryzcf
[基金申请] 国自科面上基金字体 +4 iwuli 2026-03-12 5/250 2026-03-15 17:07 by 风云无泪
[考研] 22408总分284求调剂 +3 InAspic 2026-03-13 3/150 2026-03-15 11:10 by zhq0425
[考研] 283求调剂 +6 小楼。 2026-03-12 9/450 2026-03-14 22:15 by cfxczx2015
[考研] 一志愿北京化工大学材料与化工296分求调剂 +16 稻妻小编 2026-03-09 18/900 2026-03-14 02:00 by JourneyLucky
[考研] 306求调剂 +4 唐薏薏 2026-03-09 4/200 2026-03-14 01:19 by JourneyLucky
[基金申请] 有必要更换申报口吗 20+3 fannyamoy 2026-03-11 3/150 2026-03-14 00:52 by zhanghaozhu
[考研] 279求调剂 +3 Dizzy123@ 2026-03-10 3/150 2026-03-13 23:02 by JourneyLucky
[考研] 337一志愿华南理工0805材料求调剂 +7 mysdl 2026-03-11 9/450 2026-03-13 22:43 by JourneyLucky
[考研] 材料与化工304求B区调剂 +5 邱gl 2026-03-11 6/300 2026-03-13 22:37 by JourneyLucky
[考研] 0703,333分求调剂 一志愿郑州大学-物理化学 +3 李魔女斗篷 2026-03-11 3/150 2026-03-13 22:24 by JourneyLucky
[考研] 290求调剂 +9 ADT 2026-03-11 9/450 2026-03-13 21:55 by JourneyLucky
[考研] 【0856】化学工程(085602)313 分,本科学科评估A类院校化学工程与工艺,诚求调剂 +7 小刘快快上岸 2026-03-11 7/350 2026-03-13 16:06 by ruiyingmiao
[考研] 307求调剂 +5 超级伊昂大王 2026-03-12 5/250 2026-03-13 15:56 by 棒棒球手
[考研] 328化工专硕求调剂 +4 。,。,。,。i 2026-03-12 4/200 2026-03-13 14:44 by JourneyLucky
[考研] 085600 材料与化工 295 求调剂 +10 dream…… 2026-03-10 12/600 2026-03-12 13:46 by dream……
[考研] 一志愿江南大学085701环境工程专硕总分287求调剂 +5 18266118446 2026-03-09 5/250 2026-03-11 16:51 by 2020015
[考研] 298求调剂 +3 Vv呀! 2026-03-10 3/150 2026-03-10 22:40 by 剑诗杜康
[考研] 一志愿山东大学,总分327,英语二79,有论文,有竞赛,已过四六级 +3 木木目目1 2026-03-09 3/150 2026-03-09 19:52 by yuningshan
[考研] 337求调剂 +3 睡醒,。 2026-03-09 3/150 2026-03-09 10:02 by 求调剂zz
信息提示
请填处理意见