应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜     石溪大学接受考研调剂申请>

小木虫论坛-学术科研互动平台 » 计算模拟 » beefly » 用matlab求动力学参数
【调剂】北京石油化工学院2024年16个专业接受调剂
51/1返回列表
查看: 664  |  回复: 4
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

自制奶酪

铜虫 (小有名气)

[求助] 用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
回复此楼

» 猜你喜欢
简单快乐
1楼2021-01-05 20:42:15
已阅   回复此楼   关注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的回帖

流浪在地球

新虫 (初入文坛)
您好

发自小木虫Android客户端
回复此楼
5楼2022-11-17 13:58:31
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 5 个回答
普通表情 龙 兔 虎 猫 SCI英文润色翻译-EnPapers美团队-先修改后付款
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[育儿交流] 小孩7岁，上一年级，不自信，学习较差 10+5 ljaxida 2024-04-16 26/1300 2024-04-18 13:02 by ljaxida
[教师之家] 美国教授开除中国研究生：我就不该录取你！ +5 苏东坡二世 2024-04-15 5/250 2024-04-18 08:31 by David1982
[基金申请] 申请省自然科学基金，研究区能否是省外区域 100+3 喜欢兔兔的我 2024-04-15 11/550 2024-04-17 23:49 by 喜欢兔兔的我
[考研] 275求调剂 +8 调剂0856 2024-04-14 10/500 2024-04-17 23:23 by 华北刘兵
[基金申请] 迟国泰通过向学生发放劳务费再回收的方式套取科学基金重点项目 +6 babu2015 2024-04-13 7/350 2024-04-16 20:32 by sundiv
[考研] 294求调剂 +3 694062003 2024-04-15 4/200 2024-04-16 15:01 by 邹邹哈哈
[考研] 广州大学光电信息工程专业调剂，招收物理学专业学生 +5 txhx4010 2024-04-14 7/350 2024-04-16 10:52 by domax
[考研] 312求调剂 +3 Lauhalo 2024-04-15 3/150 2024-04-16 10:16 by 19862091
[考研] 食品工程考数学283还能有什么调剂选择 +7 宇宙无敌美少女 2024-04-12 8/400 2024-04-16 10:11 by cyclodextrin726
[考研] 材料专硕329调剂遗留难民 +9 Kaylawander 2024-04-13 9/450 2024-04-15 19:21 by mumin1990
[考研] 334求调剂 +4 学药救人 2024-04-14 4/200 2024-04-15 15:05 by hunanzang
[考研] 材料与化工267求调剂 +5 工科材料267 2024-04-15 8/400 2024-04-15 14:11 by wyx颖颖
[考研] 287求调剂 +6 南沨 2024-04-14 6/300 2024-04-14 23:08 by lincunhui
[考研] 0703化学300求调剂 +19 @Xqy668800 2024-04-11 20/1000 2024-04-14 22:56 by 永字号
[考研] 一志愿大连理工大学342光学（070207），不挑学校希望老师捞一把！ +3 wxw_615 2024-04-12 4/200 2024-04-14 20:30 by zw_muchong
[考研] 334求调剂 +4 学药救人 2024-04-13 6/300 2024-04-13 20:27 by 献世的王
[考研] 339求调剂 +9 材料调剂111 2024-04-12 10/500 2024-04-13 17:07 by lincunhui
[考研] 314求调剂 +18 阿洛法 2024-04-12 20/1000 2024-04-13 13:40 by kxxmaqq
[考研] 北信科第二批调剂，工科相关 +4 biser 科研人 2024-04-11 6/300 2024-04-12 21:26 by biser 科研人
[考研] 267求调剂 +5 Wang可 2024-04-11 5/250 2024-04-12 19:29 by changhezrp
信息提示
关闭
请填处理意见
关闭