 |
 |
 |
 |
 |
 |
||
 |
 |
 |
 |
 |
小木蟲論壇-學術科研互動平臺 » 計算模擬區(qū) » 仿真模擬 » MATLAB » 【求助】修改程序
| 3 | 1/1 | 返回列表 |
| 查看: 590 | 回復: 2 | |||
[交流]
【求助】修改程序
|
|
目前下面的程序能得到r=kPA^ n PB^ m中的k、n、m值，但是我想把程序修改成能直接求r=k0e(-Ea/RT)PA^nPB^m，中的k0、Ea、n和m。應該如何修改程序，求大神賜教，不甚感激。 數據如下： T=310+273.15K時 PA0 = [3.7449 3.7500 3.7548 3.7568]; PB0 = [0.0807 0.0835 0.0861 0.0871]; r = [0.0019 0.0022 0.0024 0.0025]; T=330+273.15K時 PA0 = [3.7138 3.7215 3.7448 3.7572]; PB0 = [0.04 0.0458 0.0618 0.071]; r = [0.0041 0.0097 0.0109 0.0111]; Ps：R是常數R=8.3145 r=kPA^nPB^m的程序 function KineticsEst3 % 動力學參數辨識: 用微分法進行反應速率分析得到速率常數k和反應級數n % Reaction of the type -- rate = kCA^order % order - reaction order % rate -- reaction rate vector % CA -- concentration vector for reactant A % T -- vector of reaction time % N -- number of data points % k- reacion rate constant clear all clc global negr0m PA0 PB0 PA0 = [3.7449 3.7500 3.7548 3.7568]; PB0 = [0.0807 0.0835 0.0861 0.0871]; % negr0m = [0.5 0.63 0.83 1.0 1.28 0.33 0.8 1.5 2.21 3.33]; % Problems ... negr0m = [0.0019 0.0022 0.0024 0.0025]; % 用多變量線性回歸方法估計動力學參數 y = log(negr0m); x1 = log(PA0); x2 = log(PB0); y = y'; X = [ones(size(y)) x1' x2']; [b,bint] = regress(y,X,0.1); k = exp(b(1)) n = b(2) m = b(3) % 用多變量非線性回歸方法（以線性回歸的結果作為非線性回歸的初值） beta0 = [k n m]; lb = [0 0 0]; ub = [1 3 3]; [beta,resnorm,resid,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,beta0,lb,ub); ci = nlparci(beta,resid,jacobian) % 殘差關于擬合值的殘差圖 negrA0c = Rate(beta,PA0,PB0); plot(negrA0c,resid,'*') xlabel('反應速率擬合值, torr s^-^1') ylabel('殘差R, torr s^-^1') refline(0,0) % 參數辨識結果 fprintf('Estimated Parameters:\n') fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1)) fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2)) fprintf('\tm = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3)) fprintf('\tThe sum of the squares is: %.1e\n\n',resnorm) % ------------------------------------------------------------------ function f = ObjFunc(beta) global negr0m PA0 PB0 f = negr0m - Rate(beta,PA0,PB0); % ------------------------------------------------------------------ function negrA = Rate(beta,PA,PB) negrA = beta(1)*PA.^beta(2).*PB.^beta(3); % k=beta(1),n=beta(2),m=beta(3); |
回復此樓» 猜你喜歡
08工學調劑
已經有3人回復
-
一志愿吉林大學材料學碩321求調劑
已經有13人回復
-
一志愿西南交通 專碩 材料355 本科雙非 求調劑
已經有3人回復
-
295復試調劑
已經有4人回復
-
298-一志愿中國農業(yè)大學-求調劑
已經有8人回復
-
材料學碩297已過四六級求調劑推薦
已經有4人回復
-
能源材料化學課題組招收碩士研究生8-10名
已經有14人回復
-
0856調劑，是學校就去
已經有8人回復
-
材料080500調劑求收留
已經有5人回復
-
271材料工程求調劑
已經有7人回復
» 搶金幣啦！回帖就可以得到:
-
電力全國重點實驗室雙一流A類長江學者團隊招2026年全日制博士1名（4月3日截止）
+2/100
-
招收材料與化學相關專業(yè)碩士研究生
+1/96
-
濟南大學林秀娟教授招收材料類博士生
+1/86
-
春天就該逛吃逛吃
+1/85
-
中國藥科大學211-課題組招收2026博士/碩士
+1/53
-
華南師范大學（211）博士招生- 電子、自動化、機械、生物學、物理相關專業(yè)
+2/34
-
招聘科研助理——儲能方向（儲能、流體傳熱、CFD模擬）
+1/26
-
招聘科研助理——儲能方向（儲能、流體傳熱、CFD模擬）
+1/25
-
26年博士自薦
+1/16
-
重慶文理學院化學與環(huán)境工程學院 2026 屆碩士研究生調劑通知
+1/15
-
北京航空航天大學招生微流控與智能打印技術方向博士研究生
+2/14
-
南京醫(yī)科大學-膜蛋白結構與功能課題組-招收2026級博士研究生
+1/10
-
新加坡國立大學藥學系化學生物學課題組招PhD
+1/8
-
【博士招生】武漢科技大學招收材料、化工、環(huán)境、冶金類2026年“申請-考核”制博士生
+1/7
-
招收調劑生 國家重點實驗室（廣西某師范大學、B區(qū)）
+1/7
-
多倫多大學 & Synetic Labs 招博士后/研究員
+1/6
-
中國科學院過程工程研究所助理研究員（正式職工）和博士后招聘啟事
+1/4
-
招收化工與材料學科點研究生
+1/4
-
“機械、材料”招生
+1/3
-
中山大學醫(yī)學院趙博課題組碩士招生
+1/1
★ ★
小扶搖(金幣+10):太謝謝了，我先研究下，有問題的話還要麻煩下 2010-12-11 09:38:01
zzuwangshilei(金幣+2):積極參與,辛苦了 2010-12-11 10:36:53
小扶搖(金幣+10):太謝謝了，我先研究下，有問題的話還要麻煩下 2010-12-11 09:38:01
zzuwangshilei(金幣+2):積極參與,辛苦了 2010-12-11 10:36:53
|
我把你的程序修改了一下，能運行了，但是需要增加一組實驗數據，如果不懂，請給我短消息~ function KineticsEst3_1 % 動力學參數辨識: 用微分法進行反應速率分析得到速率常數k和反應級數n % Reaction of the type -- rate = kCA^order % order - reaction order % rate -- reaction rate vector % CA -- concentration vector for reactant A % T -- vector of reaction time % N -- number of data points % k- reacion rate constant clear all clc global negr0m PA0 PB0 RT PA0 = [3.7449 3.7500 3.7548 3.7568 3.7588 ]; PB0 = [0.0807 0.0835 0.0861 0.0871 0.0891]; RT=8.3*(310+273)*ones(1,5); % negr0m = [0.5 0.63 0.83 1.0 1.28 0.33 0.8 1.5 2.21 3.33]; % Problems ... negr0m = [0.0019 0.0022 0.0024 0.0025 0.0028]; % 用多變量線性回歸方法估計動力學參數 y = log(negr0m); x1 = log(PA0); x2 = log(PB0); x3 =log(8.3*(310+273))*ones(1,5); y = y'; X = [ones(size(y)) x1' x2' x3']; [b,bint] = regress(y,X,0.1); k0 = exp(b(1)) n = b(2) m = b(3) Ea = b(4) % 用多變量非線性回歸方法（以線性回歸的結果作為非線性回歸的初值） beta0 = [k0 n m Ea]; lb = [0 0 0 -10]; ub = [1 3 3 10]; [beta,resnorm,resid,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,beta0,lb,ub); ci = nlparci(beta,resid,jacobian) % 殘差關于擬合值的殘差圖 negrA0c = Rate(beta,PA0,PB0,RT); plot(negrA0c,resid,'*') xlabel('反應速率擬合值, torr s^-^1') ylabel('殘差R, torr s^-^1') refline(0,0) % 參數辨識結果 fprintf('Estimated Parameters:\n') fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1)) fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2)) fprintf('\tm = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3)) fprintf('\tEa = %.2f ± %.2f\n',beta(4),ci(4,2)-beta(4)) fprintf('\tThe sum of the squares is: %.1e\n\n',resnorm) %------------------------------------------------------------------ function f = ObjFunc(beta) global negr0m PA0 PB0 RT f = negr0m - Rate(beta,PA0,PB0,RT); % ------------------------------------------------------------------ function negrA = Rate(beta,PA,PB,RT) negrA = beta(1)*PA.^beta(2).*PB.^beta(3).*exp(beta(4)./RT); % k=beta(1),n=beta(2),m=beta(3);3.75883.7588 |
2樓2010-12-11 09:28:21
robert2020:建議使用“引用回復該帖”，不然對方收不到你的信息。 2010-12-16 09:15:35
|
clear all clc global negr0m PA0 PB0 RT PA0 = [3.7377 3.7449 3.7500 3.7548 3.7568 3.5690 3.5859 3.5970 3.6046 3.6119 3.6161]; PB0 = [0.0769 0.0807 0.0835 0.0861 0.0871 0.0389 0.0468 0.0520 0.0554 0.0581 0.0601]; RT1=8.3145*(310+273.15)*ones(1,5); RT2=8.3145*(330+273.15)*ones(1,6); RT=[RT1 RT2]; % negr0m = [0.5 0.63 0.83 1.0 1.28 0.33 0.8 1.5 2.21 3.33]; % Problems ... negr0m = [0.0023 0.0024 0.0024 0.0025 0.0025 0.0040 0.0042 0.0044 0.0046 0.0048 0.0049]; % 用多變量線性回歸方法估計動力學參數 y = log(negr0m); x1 = log(PA0); x2 = log(PB0); x3 =log(RT); y = y'; X = [ones(size(y)) x1' x2' x3']; [b,bint] = regress(y,X,0.1); k0 = exp(b(1)) n = b(2) m = b(3) Ea = b(4) % 用多變量非線性回歸方法（以線性回歸的結果作為非線性回歸的初值） beta0 = [k0 n m Ea]; lb = [0 0 0 -10]; ub = [1 3 3 10]; [beta,resnorm,resid,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,beta0,lb,ub); ci = nlparci(beta,resid,jacobian) % 殘差關于擬合值的殘差圖 negrA0c = Rate(beta,PA0,PB0,RT); plot(negrA0c,resid,'*') xlabel('反應速率擬合值, torr s^-^1') ylabel('殘差R, torr s^-^1') refline(0,0) % 參數辨識結果 fprintf('Estimated Parameters:\n') fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1)) fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2)) fprintf('\tm = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3)) fprintf('\tEa = %.2f ± %.2f\n',beta(4),ci(4,2)-beta(4)) fprintf('\tThe sum of the squares is: %.1e\n\n',resnorm) %------------------------------------------------------------------ function f = ObjFunc(beta) global negr0m PA0 PB0 RT f = negr0m - Rate(beta,PA0,PB0,RT); % ------------------------------------------------------------------ function negrA = Rate(beta,PA,PB,RT) negrA = beta(1)*PA.^beta(2).*PB.^beta(3).*exp(-beta(4)./RT); % k=beta(1),n=beta(2),m=beta(3);Ea=beta(4); 為什么最終的結果和 lb = [0 0 0 -10]; ub = [1 3 3 10]; 有很大關系，該如何解決 [ Last edited by 小扶搖 on 2010-12-15 at 13:35 ] |
3樓2010-12-15 13:30:40
相關版塊跳轉
我要訂閱樓主 小扶搖 的主題更新
| 3 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 一志愿吉林大學材料學碩321求調劑 +10 | Ymlll 2026-03-18 | 13/650 |
|
|---|---|---|---|---|
|
[考研] 298-一志愿中國農業(yè)大學-求調劑 +8 | 手機用戶 2026-03-17 | 8/400 |
|
|
|
[考研] 一志愿武漢理工材料工程專碩調劑 +6 | Doleres 2026-03-19 | 6/300 |
|
|
|
[考研] 304求調劑 +5 | 曼殊2266 2026-03-18 | 5/250 |
|
|
|
[考研] 295材料求調劑，一志愿武漢理工085601專碩 +3 | Charlieyq 2026-03-19 | 3/150 |
|
|
|
[考研] 復試調劑 +4 | z1z2z3879 2026-03-14 | 6/300 |
|
|
|
[考研] 一志愿西安交通大學材料工程專業(yè) 282分求調劑 +5 | 楓橋ZL 2026-03-18 | 7/350 |
|
|
|
[考研] 求調劑 +3 | Mqqqqqq 2026-03-19 | 3/150 |
|
|
|
[考研] 287求調劑 +3 | 晨昏線與星海 2026-03-19 | 4/200 |
|
|
|
[教師之家] 焦慮 +9 | 水冰月月野兔 2026-03-13 | 13/650 |
|
|
|
[考研] 材料專碩306英一數二 +10 | z1z2z3879 2026-03-16 | 13/650 |
|
|
|
[考研] 301求調劑 +9 | yy要上岸呀 2026-03-17 | 9/450 |
|
|
|
[考研] 326求調劑 +5 | 上岸的小葡 2026-03-15 | 6/300 |
|
|
|
[考研] 275求調劑 +4 | 太陽花天天開心 2026-03-16 | 4/200 |
|
|
|
[考研] 藥學383 求調劑 +3 | 藥學chy 2026-03-15 | 4/200 |
|
|
|
[考研] 機械專碩325，尋找調劑院校 +3 | y9999 2026-03-15 | 5/250 |
|
|
|
[考研] 求老師收留調劑 +4 | jiang姜66 2026-03-14 | 5/250 |
|
|
|
[考研] 070305求調劑 +3 | mlpqaz03 2026-03-14 | 4/200 |
|
|
|
[考研] 288求調劑 +4 | 奇點0314 2026-03-14 | 4/200 |
|
|
|
[考研] 復試調劑 +3 | 呼呼？~+123456 2026-03-14 | 3/150 |
|