版塊導(dǎo)航: 正在加載中...

應(yīng)《網(wǎng)絡(luò)安全法》要求，自2017年10月1日起，未進(jìn)行實(shí)名認(rèn)證將不得使用互聯(lián)網(wǎng)跟帖服務(wù)。為保障您的帳號(hào)能夠正常使用，請(qǐng)盡快對(duì)帳號(hào)進(jìn)行手機(jī)號(hào)驗(yàn)證，感謝您的理解與支持！

24小時(shí)熱門版塊排行榜

返回列表

_romantic_鎮(zhèn)

木蟲 (著名寫手)

應(yīng)助: 2 (幼兒園)
金幣: 3339
散金: 387
紅花: 6
帖子: 1221
在線: 677小時(shí)
蟲號(hào): 3186948
注冊(cè): 2014-05-07
專業(yè): 化學(xué)工程及工業(yè)化學(xué)

[求助] 初學(xué)者求助，matlab模仿程序已有1人參與

function KineticsEst5
% 動(dòng)力學(xué)ODE方程模型的參數(shù)估計(jì)

clear all
clc

k0 = [0  0  0  0];       % 參數(shù)初值
lb = [ -inf  -inf  -inf  -inf ];                % 參數(shù)下限
ub = [+inf  +inf  +inf  +inf ]; % 參數(shù)上限
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: 實(shí)驗(yàn)數(shù)據(jù)[x1 x2 x3]

% 使用函數(shù)fmincon()進(jìn)行參數(shù)估計(jì)
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函數(shù)fmincon()估計(jì)得到的參數(shù)值為:\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;

% 使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì)
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\n')
Output

% 以函數(shù)fmincon()估計(jì)得到的結(jié)果為初值，使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì)
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()的結(jié)果為初值，使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\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) ;

% 以函數(shù)fmincon()估計(jì)得到的結(jié)果為初值，使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì)
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()的結(jié)果為初值，使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\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))
];

回復(fù)此樓

» 猜你喜歡

西安電子科學(xué)大學(xué)杭州研究院劉麗香教授招收智能多模態(tài)傳感器和微型儲(chǔ)能器件方向博士已經(jīng)有13人回復(fù)
重慶交通大學(xué)光子學(xué)微結(jié)構(gòu)與器件課題組招收2026年碩士研究生信息已經(jīng)有1人回復(fù)
物理學(xué)I論文潤(rùn)色/翻譯怎么收費(fèi)? 已經(jīng)有204人回復(fù)
一志愿鄭大材料學(xué)碩298分，求調(diào)劑已經(jīng)有6人回復(fù)
尋合作：應(yīng)力腐蝕多尺度模擬已經(jīng)有3人回復(fù)
考研交流已經(jīng)有0人回復(fù)
【新加坡】納米電子器件項(xiàng)目組有“聯(lián)合培養(yǎng)博士生”名額已經(jīng)有0人回復(fù)
【新加坡】納米電子器件項(xiàng)目組有“聯(lián)合培養(yǎng)博士生”名額已經(jīng)有0人回復(fù)

1樓 2021-12-20 20:54:40

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

hzlhm

至尊木蟲 (著名寫手)

應(yīng)助: 387 (碩士)
金幣: 17807
紅花: 53
帖子: 2879
在線: 606.3小時(shí)
蟲號(hào): 1517335
注冊(cè): 2011-11-30
性別: GG
專業(yè): 常微分方程與動(dòng)力系統(tǒng)

【答案】應(yīng)助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
_romantic_鎮(zhèn): 金幣+50, ★有幫助 2022-03-15 11:07:58

程序中的問題：
1、tspan = [0.00 : 0.01 : 0.44]與tspan=[0 0.222 0.333 0.444]沒有對(duì)應(yīng)
2、f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2) 后面的代碼是多余的，計(jì)算的結(jié)果不使用
執(zhí)行結(jié)果：

使用函數(shù)fmincon()估計(jì)得到的參數(shù)值為:
k1 = -0.6403
k2 = 1.8125
k3 = 3.6968
k4 = 1.9950
The sum of the squares is: 4.0

贊一下(1人)

回復(fù)此樓

QQ：2120156492

2樓2021-12-27 22:26:21

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

相關(guān)版塊跳轉(zhuǎn) 我要訂閱樓主 _romantic_鎮(zhèn) 的主題更新