查看: 340  |  回复: 3

自制奶酪

铜虫 (小有名气)

[求助] 用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的回帖

独孤神宇

版主 (职业作家)

【答案】应助回帖

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

lsqnonlin 函数 ObjFunc4LNL 缺失
数值计算
2楼2021-01-05 21:03:26
已阅   回复此楼   关注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的回帖

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的回帖
相关版块跳转 我要订阅楼主 自制奶酪 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 建议JJW建立“行为不端举报奖励机制”,包括基金倾斜 +9 hneyxiong 2022-08-16 10/500 2022-08-16 12:36 by kinlin13
[基金申请] 啥时候公示? +8 nickzhaisdnu 2022-08-15 10/500 2022-08-16 12:12 by sansan99
[基金申请] 投票:  2022年8月10日,申请书的文件名来投票啦~~~~~~ +24 mmffdd26 2022-08-10 28/1400 2022-08-16 12:08 by derickofaze
[基金申请] 今天我生日,青基祈福 +38 plumage5 2022-08-16 38/1900 2022-08-16 12:00 by B_B
[基金申请] 我这个研究基础申请省自然科学基金面上项目有戏吗? +12 yue_s008 2022-08-15 14/700 2022-08-16 11:48 by yue_s008
[硕博家园] 社恐怎么和老师正常交流 +12 学术菜菜鸟 2022-08-15 23/1150 2022-08-16 09:32 by 学术菜菜鸟
[基金申请] 请问现在可以用国际合作项目查结果吗? +9 猪猪迎 2022-08-16 10/500 2022-08-16 09:25 by wmfsnow
[基金申请] 第四年申请青基了,希望有好结果 +3 一川烟草 2022-08-16 3/150 2022-08-16 08:55 by flybigcat
[基金申请] pdf变了 +28 无情的雨2011 2022-08-15 29/1450 2022-08-16 07:39 by 小石头0125
[文学芳草园] 平静 +3 myrtle 2022-08-10 6/300 2022-08-16 06:35 by 飘过的晚辈
[基金申请] 我们知道科研之友阅读量增加与中不中没关系,为什么评审阅读量会明显增加呢 +6 翰海2022 2022-08-15 14/700 2022-08-15 23:19 by ryzzwmc
[基金申请] 现在申请书名称都是十位数字-proposal吗? +10 jane_012 2022-08-14 12/600 2022-08-15 20:56 by jane_012
[基金申请] 基金编号加proposal是没中吗? +11 fengjinglove 2022-08-15 13/650 2022-08-15 17:46 by cfkybfq
[找工作] 还有不坑的学校吗? +13 菜鸟夭夭 2022-08-12 13/650 2022-08-15 00:42 by shenqi2
[论文投稿] 投稿当天就Required reviews completed +5 wfyuan1975 2022-08-12 8/400 2022-08-14 22:45 by kingvsqueenx
[论文投稿] 对自己科研能力不自信,结果出现奇迹 +8 Xzhang-23 2022-08-10 10/500 2022-08-12 20:04 by 选好路
[公派出国] +6 taotaol 2022-08-11 12/600 2022-08-12 14:14 by misterfeng
[基金申请] 文件后缀变化了626efXXXXX-ProposalPDF-202202 +13 H098 2022-08-11 16/800 2022-08-11 15:36 by libeny
[基金申请] 那些花儿! +6 了却无痕 2022-08-10 16/800 2022-08-11 09:17 by 医学生的成长
[教师之家] 南京晓庄 金陵科技待遇 +4 jamie1115 2022-08-10 5/250 2022-08-10 17:46 by owl2010
信息提示
请填处理意见