版塊導(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í)熱門(mén)版塊排行榜

返回列表

805592836_01

鐵蟲(chóng) (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 823
帖子: 10
在線: 23.4小時(shí)
蟲(chóng)號(hào): 2393210
注冊(cè): 2013-04-01
性別: MM
專(zhuān)業(yè): 能源化工

[求助] 未完成的matlab 已有2人參與

有實(shí)驗(yàn)測(cè)得的一系列時(shí)間點(diǎn)
t=（0,180,300,420,600,900,1200,1800,2400）
CA=（1.845，1.414,1.237,1.065,0.873,0.629,0.440,0.226,0.124）

微分方程：k(1)*k(2)*k(3)*CA0*(CB.^3)./(1+k(2)*CA+k(3)^(1/6)*CB.^(1/2)).^7
CB=0.4為定值
希望由此得到三個(gè)參數(shù)k(1),k(2),k(3)，以及實(shí)驗(yàn)值與擬合值的比較圖。現(xiàn)已經(jīng)編了部分代碼，但不是很確定對(duì)否，希望高手可以修改一下

clear all; clc
k0=[1,1,1];
lb=[0,0,0];
ub=[+inf,+inf,+inf];

tspan=[0,180,300,420,600,900,1200,1800,2400];

yexp=[1.414,1.237,1.065,0.873,0.629,0.440,0.226,0.124]';
y0=1.845;

% 使用fmincon()進(jìn)行參數(shù)估計(jì)
[k,fval,flag]=fmincon(@ObjFunc4Fmincon,k0,lb,ub,option,tspan,y0,yexp);
fprintf('\tk1=%.4f\n',k(1)),  %本征動(dòng)力學(xué)參數(shù)
fprintf('\tk2=%.4f\n',k(2)),  % A的吸附平衡常數(shù)KA
fprintf('\tk3=%.4f\n',k(3)),  % H2的吸附平衡常數(shù)KB
fprintf('The sum of the squares is:%.1e\n\n',fval),
k_fmincon=k;

% 使用函數(shù)Isqnonlin()進(jìn)行參數(shù)估計(jì)
yy=[y0 yexp'];
options = optimset('MaxFunEvals',100000)
[k,resmorm,residual,exitflag,output,lambda,jacobian]=...
Isqnonlin(@ObjFunc4LNL,k0,lb,ub,options,tspan,y0,yexp);
ci=nlparci(k,residual,jacobian);
fprintf('n\n使用Iaqunonlin()估計(jì)得到的參數(shù)值:\n'),Output

% 以函數(shù)fmincon()估計(jì)得到的結(jié)果為初值，使用Isqnonlin()進(jìn)行參數(shù)估計(jì)
k0=k_fmincon;
[k,resmorm,residual,exitflag,output,lambda,jacobian]=...
Isqnonlin(@ObjFunc4LNL,k0,lb,ub,options,tspan,y0,yexp);
ci=nlparci(k,residual,jacobian);
fprintf('n\n以fmincon()結(jié)果為初值，使用函數(shù)Isqnonlin()估計(jì)得到的參數(shù)值為:\n')
Output

%------------------------------------------------------------

function f=ObjFunc4Fmincon(k,tspan,y0,yexp)
[t y]=ode45(@KineticsEqs,tspan,y0,yexp,k);

%----------------------------------------------------------

function f=Objfunc4LNL(k,tspan,yexp)
[t y]=ode45(@KineticEqs,tspan,tspan,[],k);

%----------------------------------------------------------
function dydt = KineticsEqs(t,y,k)
y0=1.845;
CB=y;
CB0=0.4;
CA0=y0;
k1=k(1);
KA=k(2);
KB=k(3);
dxdt=k(1)*k(2)*k(3)*CA0*(CB.^3)./(1+k(2)*CA0+k(3)^(1/6)*CB.^(1/2)).^7

回復(fù)此樓

» 猜你喜歡

哈爾濱理工大學(xué)物理系招收物理學(xué)考研調(diào)劑已經(jīng)有6人回復(fù)
0702一志愿吉大B區(qū)求調(diào)劑已經(jīng)有5人回復(fù)
物理學(xué)I論文潤(rùn)色/翻譯怎么收費(fèi)? 已經(jīng)有194人回復(fù)
求調(diào)劑已經(jīng)有0人回復(fù)
0702一志愿吉大B區(qū)求調(diào)劑有論文已經(jīng)有0人回復(fù)
請(qǐng)問(wèn)還有沒(méi)有用Latex寫(xiě)文章的小伙伴們？已經(jīng)有0人回復(fù)
光學(xué)工程學(xué)碩調(diào)劑信息已經(jīng)有31人回復(fù)
歡迎加入課題組已經(jīng)有0人回復(fù)
散金幣，求好運(yùn)，祝面上順利！已經(jīng)有65人回復(fù)
華東師范大學(xué)芯片設(shè)計(jì)徐瓏真實(shí)水平咋樣已經(jīng)有1人回復(fù)
調(diào)劑已經(jīng)有0人回復(fù)

» 本主題相關(guān)價(jià)值貼推薦，對(duì)您同樣有幫助:

matlab中的計(jì)算問(wèn)題？已經(jīng)有8人回復(fù)
MATLAB求助，關(guān)于importdata的問(wèn)題已經(jīng)有6人回復(fù)
用matlab制作的M文件保存在哪里才能直接調(diào)用？已經(jīng)有9人回復(fù)
matlab一運(yùn)行就出現(xiàn)的問(wèn)題已經(jīng)有7人回復(fù)
matlab的m文件關(guān)聯(lián)問(wèn)題已經(jīng)有17人回復(fù)
求助一個(gè)matlab關(guān)于分類(lèi)討論的問(wèn)題已經(jīng)有9人回復(fù)
matlab編程中遇到的問(wèn)題已經(jīng)有4人回復(fù)
幫忙修改Matlab中出現(xiàn)的錯(cuò)誤已經(jīng)有4人回復(fù)
matlab讀取圖像的問(wèn)題已經(jīng)有12人回復(fù)
分享matlab的一些代碼（入門(mén)貼）已經(jīng)有19人回復(fù)
求MATLAB清華版的教學(xué)視頻已經(jīng)有14人回復(fù)
關(guān)于Matlab作圖的問(wèn)題已經(jīng)有14人回復(fù)
MATLAB初學(xué)者，這樣的式子怎么就是畫(huà)不了圖，總是出錯(cuò)!求指導(dǎo)!!!!!!!!!!!!! 已經(jīng)有8人回復(fù)
matlab程序?qū)D中的數(shù)據(jù)記錄下來(lái) 已經(jīng)有4人回復(fù)
Matlab在地學(xué)中的應(yīng)用已經(jīng)有24人回復(fù)
matlab中dydt的編寫(xiě)要求已經(jīng)有12人回復(fù)
求助svm在matlab中的應(yīng)用已經(jīng)有7人回復(fù)
關(guān)于matlab 2010b 的一個(gè)小問(wèn)題已經(jīng)有7人回復(fù)
MATLAB安裝后出現(xiàn)一堆java的問(wèn)題已經(jīng)有12人回復(fù)
【求助】關(guān)于matlab的高維數(shù)組已經(jīng)有6人回復(fù)

加油加油加油……

1樓 2014-03-12 08:57:58

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

dingd

鐵桿木蟲(chóng) (職業(yè)作家)

計(jì)算強(qiáng)帖: 4
應(yīng)助: 1641 (講師)
金幣: 15037.3
散金: 101
紅花: 234
帖子: 3410
在線: 1223.7小時(shí)
蟲(chóng)號(hào): 291104
注冊(cè): 2006-10-28

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

感謝參與，應(yīng)助指數(shù) +1

k1, k2, k3要求都大于0的話效果會(huì)很差的。

贊一下

回復(fù)此樓

2樓2014-03-12 09:56:13

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

805592836_01

鐵蟲(chóng) (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 823
帖子: 10
在線: 23.4小時(shí)
蟲(chóng)號(hào): 2393210
注冊(cè): 2013-04-01
性別: MM
專(zhuān)業(yè): 能源化工

引用回帖:

2樓: Originally posted by dingd at 2014-03-12 09:56:13
k1, k2, k3要求都大于0的話效果會(huì)很差的。

嗯，都要大于0的……

贊一下

回復(fù)此樓

加油加油加油……

3樓2014-03-12 10:02:43

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

dingd

鐵桿木蟲(chóng) (職業(yè)作家)

計(jì)算強(qiáng)帖: 4
應(yīng)助: 1641 (講師)
金幣: 15037.3
散金: 101
紅花: 234
帖子: 3410
在線: 1223.7小時(shí)
蟲(chóng)號(hào): 291104
注冊(cè): 2006-10-28

引用回帖:

3樓: Originally posted by 805592836_01 at 2014-03-12 10:02:43
嗯，都要大于0的……...

不可能有好結(jié)果，檢查下模型吧。

贊一下

回復(fù)此樓

4樓2014-03-12 10:35:53

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

805592836_01

鐵蟲(chóng) (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 823
帖子: 10
在線: 23.4小時(shí)
蟲(chóng)號(hào): 2393210
注冊(cè): 2013-04-01
性別: MM
專(zhuān)業(yè): 能源化工

引用回帖:

4樓: Originally posted by dingd at 2014-03-12 10:35:53
不可能有好結(jié)果，檢查下模型吧。...

嗯，好的吧！同樣感謝

贊一下

回復(fù)此樓

加油加油加油……

5樓2014-03-12 11:31:29

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

月只藍(lán)

主管區(qū)長(zhǎng) (職業(yè)作家)

專(zhuān)家經(jīng)驗(yàn): +1059
計(jì)算強(qiáng)帖: 8
應(yīng)助: 1712 (講師)
貴賓: 8.888
金幣: 68127.7
散金: 1938
紅花: 443
沙發(fā): 4
帖子: 4373
在線: 3291.4小時(shí)
蟲(chóng)號(hào): 1122189
注冊(cè): 2010-10-14
專(zhuān)業(yè): 宇宙學(xué)
管轄: 計(jì)算模擬區(qū)

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

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感謝參與，應(yīng)助指數(shù) +1
805592836_01: 金幣+15, ★★★很有幫助, 看來(lái)還是模型不合適，三個(gè)參數(shù)要都大于0，不過(guò)還是感謝大神的幫忙。 2014-03-12 21:06:07

原代碼問(wèn)題很多，比如MATLAB中不存在Isqnonlin函數(shù)，系拼寫(xiě)錯(cuò)誤，正確的為 lsqnonlin；此外用fmincon函數(shù)算出來(lái)的結(jié)果輸入給lsqnonlin函數(shù)作初值意義不是很大，單獨(dú)調(diào)用lsqnonlin即可。
根據(jù)實(shí)驗(yàn)數(shù)據(jù)可知，數(shù)據(jù)呈單調(diào)遞減的趨勢(shì)，擬合公式是關(guān)于時(shí)間的導(dǎo)數(shù)，應(yīng)小于0，所以k1 k2 k3不應(yīng)該全部大于0，否則不可能有號(hào)的擬合效果。
在以下給出程序中，把k1的取值范圍可為0，具體代碼以及計(jì)算結(jié)果如下。

CODE:

function test_123

clear all;clc

format long

tspan=[0,180,300,420,600,900,1200,1800,2400];

yexp=[1.414,1.237,1.065,0.873,0.629,0.440,0.226,0.124]';

y0=1.845;

k0=[-1 1 1];

lb=-[inf 0 0];

ub=+[inf inf inf];        

yy=[y0 yexp'];

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...

    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,y0,yexp);

ci = nlparci(k,residual,jacobian);

fprintf('\n\n使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\n')

fprintf('\tk1=%.4f\n',k(1)),  %本征動(dòng)力學(xué)參數(shù)

fprintf('\tk2=%.4f\n',k(2)),  % A的吸附平衡常數(shù)KA

fprintf('\tk3=%.4f\n',k(3)),  % H2的吸附平衡常數(shù)KB

fprintf('  The sum of the squares is: %.1e\n\n',resnorm)

ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,y0,[],k);

[ttt XXsim] = ode45(@KineticsEqs,tspan,y0,[],k); 

y=XXsim(2:end);

xexp=yexp;

R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2);

fprintf('\n\t相關(guān)系數(shù)之平方R^2 = %.6f',R2);

plot(ts,ys,'b',tspan,yy,'or'),legend('計(jì)算值','實(shí)驗(yàn)值','Location','best'),

%---------------------------------------------------------

function f = ObjFunc(k,tspan,y0,yexp)           % 目標(biāo)函數(shù)

[t Xsim] = ode45(@KineticsEqs,tspan,y0,[],k) ;

ysim = Xsim(2:end);

f=ysim-yexp;

 %----------------------------------------------------------

function dxdt = KineticsEqs(t,y,k)

y0=1.845;

CB=y;

CB0=0.4;

CA0=y0;

k1=k(1);

KA=k(2);

KB=k(3);

dxdt=k(1)*k(2)*k(3)*CA0*(CB.^3)./(1+k(2)*CA0+k(3)^(1/6)*CB.^(1/2)).^7;

使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:
k1=-1.0834
k2=0.0598
k3=5.6626
The sum of the squares is: 7.8e-003

相關(guān)系數(shù)之平方R^2 = 0.995002>>

附圖1.jpg

贊一下(1人)

回復(fù)此樓

MATLAB、MS小問(wèn)題、普通問(wèn)題請(qǐng)發(fā)帖求助！時(shí)間精力有限，恕不接受無(wú)償私信求助。

6樓2014-03-12 13:56:30

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

相關(guān)版塊跳轉(zhuǎn) 我要訂閱樓主 805592836_01 的主題更新

返回列表

普通表情龍兔虎貓

最具人氣熱帖推薦 [查看全部]		作者	回/看	最后發(fā)表

[考研] 一志愿天津大學(xué)化學(xué)工藝專(zhuān)業(yè)（081702）315分求調(diào)劑 +12	yangfz 2026-03-17	12/600	2026-03-21 03:30 by JourneyLucky
[考研] 332求調(diào)劑 +4	ydfyh 2026-03-17	4/200	2026-03-21 02:20 by JourneyLucky
[考研] 華東師范大學(xué)-071000生物學(xué)-293分-求調(diào)劑 +3	研究生何瑤明 2026-03-18	3/150	2026-03-21 01:30 by JourneyLucky
[考研] 299求調(diào)劑 +4	shxchem 2026-03-20	5/250	2026-03-21 00:24 by 星空星月
[考研] 南昌大學(xué)材料專(zhuān)碩311分求調(diào)劑 +5	77chaselx 2026-03-20	5/250	2026-03-20 23:42 by lovewei0727
[考研] 308求調(diào)劑 +3	阿姐阿姐家啊 2026-03-18	3/150	2026-03-20 23:24 by JourneyLucky
[考研] 321求調(diào)劑 +9	何潤(rùn)采123 2026-03-18	11/550	2026-03-20 23:19 by JourneyLucky
[考研] 一志愿西北大學(xué) ，070300化學(xué)學(xué)碩，總分287，雙非一本，求調(diào)劑。 +4	晨昏線與星海 2026-03-19	4/200	2026-03-20 22:15 by JourneyLucky
[考研] 一志愿武理材料工程348求調(diào)劑 +3	￣^￣゜汗 2026-03-19	4/200	2026-03-20 21:01 by zhukairuo
[考研] 281求調(diào)劑（0805） +14	煙汐憶海 2026-03-16	25/1250	2026-03-20 15:47 by yuncha
[考研] 298-一志愿中國(guó)農(nóng)業(yè)大學(xué)-求調(diào)劑 +9	手機(jī)用戶(hù) 2026-03-17	9/450	2026-03-20 14:24 by 無(wú)懈可擊111
[考研] 求調(diào)劑 +3	暗涌afhb 2026-03-16	3/150	2026-03-20 00:28 by 河南大學(xué)校友
[考研] 生物學(xué)調(diào)劑招人！！！ +3	山海天嵐 2026-03-17	4/200	2026-03-19 21:34 by 怎么釋?xiě)?/a>
[考研] 0703化學(xué) 305求調(diào)劑 +4	FY_yy 2026-03-14	4/200	2026-03-19 05:54 by anny19840123
[考研] 材料專(zhuān)碩306英一數(shù)二 +10	z1z2z3879 2026-03-16	13/650	2026-03-18 14:20 by 007_lilei
[考研] 302求調(diào)劑 +10	呼呼呼。。。。 2026-03-17	10/500	2026-03-18 12:45 by Linda Hu
[考研] 268求調(diào)劑 +6	簡(jiǎn)單點(diǎn)0 2026-03-17	6/300	2026-03-18 09:04 by 無(wú)際的草原
[考研] 材料工程專(zhuān)碩274一志愿211求調(diào)劑 +6	薛云鵬 2026-03-15	6/300	2026-03-17 11:05 by 學(xué)員h26Tkc
[考研] 26考研一志愿中國(guó)石油大學(xué)(華東)305分求調(diào)劑 +3	嘉年新程 2026-03-15	3/150	2026-03-15 13:58 by 哈哈哈哈嘿嘿嘿
[考研] 本科南京大學(xué)一志愿川大藥學(xué)327 +3	麥田耕者 2026-03-14	3/150	2026-03-14 20:04 by 外星文明