24小时热门版块排行榜    

查看: 959  |  回复: 1

_romantic_镇

木虫 (著名写手)

[求助] 初学者求助,matlab模仿程序已有1人参与

function KineticsEst5
% 动力学ODE方程模型的参数估计

clear all
clc

k0 = [0  0  0  0];         % 参数初值
lb = [ -inf  -inf  -inf  -inf ];                   % 参数下限
ub = [+inf  +inf  +inf  +inf ];    % 参数上限
x0 = [8.5 28.8 27.6];
tspan=[0 0.222 0.333 0.444];
yexp =  [8.5000   28.8000   27.6000
    4.8000   23.2000   35.3000
    4.2000   21.6000   36.5000
    4.0000   21.2000   37.3000];                 % yexp: 实验数据[x1        x2 x3]

% 使用函数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('  The sum of the squares is: %.1e\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


% ------------------------------------------------------------------
function f = ObjFunc4Fmincon(k,x0,yexp)
tspan = [0.00 : 0.01 : 0.44];
[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) ;

% 以函数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

% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)
tspan = [0.00 : 0.01 : 0.44];
[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(2)- k(1)*x(1))
   (k(1)*x(2)+(k(3)-k(2)-k(4))*x(3))
   (k(3)*x(2)-k(4)*x(3))
   ];
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzlhm

至尊木虫 (著名写手)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
_romantic_镇: 金币+50, 有帮助 2022-03-15 11:07:58
程序中的问题:
1、tspan = [0.00 : 0.01 : 0.44]与tspan=[0 0.222 0.333 0.444]没有对应
2、f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2)  + sum((y(:,3)-yexp(:,3)).^2)  后面的代码是多余的,计算的结果不使用
执行结果:

使用函数fmincon()估计得到的参数值为:
        k1 = -0.6403
        k2 = 1.8125
        k3 = 3.6968
        k4 = 1.9950
The sum of the squares is: 4.0
QQ:2120156492
2楼2021-12-27 22:26:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 _romantic_镇 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 关于硕博连读的一些疑问? +4 Lwenter 2024-06-14 5/250 2024-06-17 20:42 by Lwenter
[考博] 这个博士要读吗 +9 Sea Breeze 2024-06-16 17/850 2024-06-17 19:50 by zeolitess
[基金申请] 有人中过人文社科类的博后特助吗? +3 outsider1986 2024-06-16 4/200 2024-06-17 19:35 by outsider1986
[论文投稿] 审稿问题:为什么荧光激发波长和紫外吸收波长差的大? 10+5 sdawege 2024-06-14 10/500 2024-06-17 18:54 by HH-探针
[基金申请] 我标书代码变了 +44 学员NHuqdk 2024-06-16 58/2900 2024-06-17 17:39 by wenfengw83
[基金申请] 博后面上今天有bug可以看到是否资助? +22 lyfbangong 2024-06-12 33/1650 2024-06-17 16:17 by carolloo
[博后之家] 山东大学(青岛)“天然药物生物智造”课题组 招聘“博士后”(年薪20.4-55.6万元) +5 第二种态度 2024-06-11 11/550 2024-06-17 10:08 by 小懂事k
[找工作] 应届大龄土博教职选择,石河子大学和广西大学 10+4 茅山老道士 2024-06-10 34/1700 2024-06-16 22:05 by yugege2009
[基金申请] 博后基金,博管会会提前知道消息吗? +4 yuyiang 2024-06-13 4/200 2024-06-16 11:40 by yangyuzhong4
[基金申请] BO4的YQ答辩通知发布了吗? +6 博学笃行 2024-06-11 6/300 2024-06-15 16:04 by 悲催科研狗
[基金申请] 为什么我的博后基金还在流动站审核中?不会是学院给我卡了吧? +14 王凯12 2024-06-13 26/1300 2024-06-15 15:22 by 好人与坏人
[基金申请] 博后基金,以往的结果点不开,怎么回事呢?最后一次机会了,两次都没中前面。 +7 kyukitu 2024-06-14 13/650 2024-06-15 06:46 by 我是王小帅
[基金申请] 博士后基金需要结题吗? +8 zhouchuck 2024-06-13 8/400 2024-06-14 17:27 by liuyupu132
[考博] 申博找导师 +4 疏影横斜水清浅3 2024-06-13 6/300 2024-06-14 14:31 by zxl_1105
[基金申请] 75批博后基金 +10 kyukitu 2024-06-13 13/650 2024-06-14 10:31 by kyukitu
[论文投稿] 最近写了一篇控制优化领域的文章,可以投哪里啊?有没有水一些的期刊推荐 +7 香瓜木香 2024-06-12 13/650 2024-06-14 07:05 by 香瓜木香
[论文投稿] 文章proof要求使用机构的邮箱 5+3 不可不信缘 2024-06-11 11/550 2024-06-14 07:00 by 3001160025
[考博] 博导选择 +3 bing85977 2024-06-12 3/150 2024-06-13 15:34 by 我是邱尧
[论文投稿] 摩擦磨损论文投稿 +3 jmysan 2024-06-12 3/150 2024-06-13 08:36 by 莱茵润色
[基金申请] 博后特助这周出结果吗?往年都是啥时候啊? +13 jsqy 2024-06-12 17/850 2024-06-12 19:55 by Lynn212
信息提示
请填处理意见