24小时热门版块排行榜    

查看: 949  |  回复: 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_镇 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 24年还有机会申博吗,求博导收留 +3 鄙人康康 2024-05-23 6/300 2024-05-27 00:29 by 鄙人康康
[硕博家园] 文科博在木虫上存在感好低呀 +6 hahamyid 2024-05-25 9/450 2024-05-27 00:09 by 飘过的晚辈
[硕博家园] 要不要读博 +8 王乔木 2024-05-24 8/400 2024-05-26 23:48 by 余东东锵锵
[硕博家园] 研0 +6 控制调剂yl 2024-05-25 11/550 2024-05-26 23:37 by sakuraai
[教师之家] 被惯着的学生终究要吃大亏 +21 535743368 2024-05-24 22/1100 2024-05-26 22:56 by mmxf_2011
[考博] 国家高层次人才李兴淑教授课题组急招博士研究生 +3 宁静de夏 2024-05-20 5/250 2024-05-26 19:46 by 宁多缺毋滥
[电化学] 如何证明电极上镀上金了? 10+4 刻印时光 2024-05-22 5/250 2024-05-26 19:03 by 刻印时光
[硕博家园] 博士复试,申请成绩复核,有机会翻盘吗? +22 长海二声笑 2024-05-21 29/1450 2024-05-26 14:48 by 鱼翔浅底1
[基金申请] 科研之友阅读量近一周增加了200多。 +11 hdzw9071 2024-05-24 12/600 2024-05-26 09:38 by wanghuawei
[硕博家园] 人生 +15 暮色恋伊人 2024-05-22 15/750 2024-05-26 08:23 by elainzai
[考博] 25年博士申请 +5 zjc晨 2024-05-24 8/400 2024-05-25 22:45 by zjc晨
[教师之家] 经常觉得挺累的 +13 zylfront 2024-05-22 19/950 2024-05-25 21:08 by hjc404
[基金申请] 河北省基金 50+3 晓晓爱翠翠 2024-05-23 19/950 2024-05-25 15:34 by 972937946
[论文投稿] 为什么有的影响因子高的期刊分区不高呢? +9 安处一室 2024-05-21 9/450 2024-05-25 15:24 by dhsjg
[考博] 24级求博导 +3 Hddd9 2024-05-23 3/150 2024-05-24 18:44 by 安塔瓦拉多
[基金申请] 听说面青地E09已经送了么? +6 叉烧吃叉烧 2024-05-21 9/450 2024-05-23 12:24 by 叉烧吃叉烧
[论文投稿] Neurocomputing 外审结束 +5 mollyzhang_2003 2024-05-23 5/250 2024-05-23 12:03 by nono2009
[论文投稿] 关于通讯作者 5+4 irikiar 2024-05-21 4/200 2024-05-23 09:43 by moyoushang
[教师之家] 有没有在职教师同时做博后的? +6 克雷斯 2024-05-20 8/400 2024-05-23 08:08 by 克雷斯
[基金申请] 太诡异了,五月底还有没有送审的。。 +12 hdzw9071 2024-05-21 12/600 2024-05-21 12:43 by dxcharlary
信息提示
请填处理意见