最自私的女人

新蟲 (初入文壇)

應助: 0 (幼兒園)
金幣: 136.5
紅花: 1
帖子: 24
在線: 7小時
蟲號: 2922685
注冊: 2014-01-09
專業(yè): 化學反應工程

[求助] 動力學方程組參數求解問題已有3人參與

本人動力學方程如下：
dx1/dt=k21x2+k31x3+k41x4+k51x5+k61x6+k71x7+k81x8
dx2/dt=-k21x2+k32x3+k42x4+k52x5+k62x6+k72x7+k82x8
dx3/dt=-k31x3-k32x3+k43x4+ k53x5+k63x6+k73x7+k83x8
dx4/dt=-k41x4-k42x4-k43x4+ k54x5+k64x6+k74x7+k84x8
dx5/dt=-k51x5-k52x5-k53x5-k54x5-k56x52-k57x52+k65x6+k75x7+k85x8
dx6/dt=-k61x6-k62x6-k63x6-k64x6-k65x6-k67x62+k56x52+k76x7+k86x8
dx7/dt=-k71x7-k72x7-k73x7-k74x7-k75x7-k76x7-k78x72+k57x52+k67x62+k87x8
dx8/dt=-k81x8-k82x8-k83x8-k84x8-k85x8-k86x8-k87x8-k89x82+k78x72
dx9/dt=k89x82
方程中有多個縮合反應為表觀二級（即平方），已上傳圖片供參考
k為動力學參數，k大于0，共33個
x為質量分數，x＞0，九個x之和為1
已知數據如下：
溫度時間 x1 x2 x3 x4 x5 x6 x7 x8 x9
430 0.00 0.0192 0.0236 0.1006 0.3284 0.1647 0.1327 0.1981 0.0319 0.0008
430 20.00 0.0287 0.0503 0.1313 0.3156 0.1391 0.1206 0.1755 0.0366 0.0022
430 40.00 0.0349 0.0858 0.1827 0.3058 0.1203 0.0805 0.1476 0.0384 0.0039
430 60.00 0.0423 0.1250 0.2086 0.2950 0.1005 0.0624 0.1180 0.0428 0.0053
本人想編寫一個MATLAB程序實現(xiàn)以下功能：
設k初值均為0，求解方程中的33個參數，將求得的參數代入方程中，根據t=0min的數據求得t=20、40、60min的x計算值，將x計算值與x實驗值比較，如果誤差＜5%，得到的參數就是最終參數，如果誤差大于5%，就將得到的參數最為初值循環(huán)上述計算，直到誤差＜5%。
哪位大俠能幫忙編寫一下，小女子不勝感激呢~

動力學方程組.jpg

回復此樓

» 猜你喜歡

2026年循環(huán)經濟功能材料國際會議（ICFMCE 2026）已經有0人回復
2026年第五屆電氣、電子與信息工程國際會議（ISEEIE 2026）已經有0人回復
物理學I論文潤色/翻譯怎么收費? 已經有111人回復
哈爾濱理工大學物理系招收物理學考研調劑已經有5人回復
0702一志愿吉大B區(qū)求調劑已經有5人回復
求調劑已經有0人回復
0702一志愿吉大B區(qū)求調劑有論文已經有0人回復
請問還有沒有用Latex寫文章的小伙伴們？已經有0人回復
光學工程學碩調劑信息已經有29人回復
歡迎加入課題組已經有0人回復
散金幣，求好運，祝面上順利！已經有52人回復

» 本主題相關價值貼推薦，對您同樣有幫助:

使用matlab最優(yōu)化方法擬合獲得多個動力學參數中的問題已經有4人回復
有關偽一級動力學和偽二級動力學方程已經有9人回復
用matlab求解方程已經有6人回復
求助動力學方程？？？已經有3人回復
求助用matlab擬合動力學方程已經有9人回復
如何求解帶有固定參數的超定方程組？已經有21人回復
求助，matlab自定義微分方程擬合實驗數據來求方程中的參數已經有4人回復
微分方程組求參數問題，求高人指教，感謝已經有11人回復
微分方程組參數擬合的問題（多參數動力系統(tǒng)參數估計）已經有11人回復
【求助】非線性方程組的求解問題已經有6人回復
求助matlab擬合方程組已經有12人回復
動力學方程參數估計方法已經有14人回復
SCR反應動力學方程求解的問題已經有4人回復
動力學參數擬合已經有26人回復
請教matlab反應動力學參數估計遇到的問題，謝謝已經有15人回復
車輛動力學模型的解法已經有14人回復
matlab m文件運行出錯已經有4人回復
matlab-常微分方程參數估計已經有12人回復
lsqnonlin函數擬合微分方程組參數擬合問題已經有10人回復
求回歸方程參數估計的最優(yōu)化算法matlab代碼已經有9人回復
matlab非線性參數擬合問題已經有7人回復
【求助】matlab 求解微分方程中的未知參數已經有20人回復
【求助】多元方程組求解及擬合已經有9人回復
【求助】催化反應動力學已經有5人回復
【求助】使用Matlab預估動力學方程問題已經有13人回復

1樓 2014-09-02 13:55:16

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

最自私的女人

新蟲 (初入文壇)

應助: 0 (幼兒園)
金幣: 136.5
紅花: 1
帖子: 24
在線: 7小時
蟲號: 2922685
注冊: 2014-01-09
專業(yè): 化學反應工程

引用回帖:

2樓: Originally posted by 月只藍 at 2014-09-02 14:26:53
參數估計問題或者叫常微分方程組參數擬合。參數太多，MATLAB求解估計困難不小，而且數據僅有4組。
建議你找有高版本1stOpt軟件的專家?guī)湍阕鲆幌�，推薦本版塊的專家dingd。

你好，能麻煩你幫我看一下我的下面這個程序為什么運行不了嘛？
function fourLumpK_VR440
clc;
clear all;
format long;
fourLumpData=[...
0 69.33 21.33 7.66 1.76
15 52.96 26.45 17.84 3.18
30 49.9 27.61 18.29 3.92
45 39.03 28.51 27.93 5
60 33.45 25.4 34.49 7.1
90 25.76 24.99 42.38 7.48
];
% vector of yields/%
Xexp=fourLumpData(2:6,2:5);
% oil residence time,unit: min
time= fourLumpData(2:6,1)';

% initiate parameter for optimization
X0=fourLumpData(1,2:5);
% k1 k2 k3 k4 k5
K0 =[0.0014 0.0012 0.0005 0.0020 0.0013];
lb=0;
ub=inf;

%計時開始擬合
tic;

K=lsqnonlin(@objFun,K0,lb,ub,[],X0,Xexp,time);
% 擬合結果顯示
disp(K)
% 計時結束
disp('耗時：');
toc;
% 保存擬合結果
save fourLumpK_VR440;

% 檢驗-----------------------------------------------
tf=max(time);
tspan=[0:1:tf];
[nr,nc]=size(Xexp);
XcalFinal=zeros(nr,nc);
for i=1:nr
[t,Xcal]=ode45(@modelEquation,tspan,X0,[],K,tf);
n=find(t>=time(i));
m=n(1);
%線性內插
Xc=Xcal(m,

-(Xcal(m,

-Xcal(m-1,

)/(t(m)-t(m-1))*(t(m)-time(i));
XcalFinal(i,

=Xc;
save XcalFinal
end
objX = abs((XcalFinal - Xexp)./100).^2;
F=sum(sum(objX))^0.5;
error=(XcalFinal - Xexp)
avrError=sum(abs(error))/nr
relatError=(XcalFinal - Xexp)./Xexp*100
avrRelError=sum(abs(relatError))/nr

rou_2=1-sum((XcalFinal - Xexp).^2)./sum(Xexp.^2)
Fc_11=(sum(XcalFinal.^2)-sum((XcalFinal - Xexp).^2))/5./(sum((XcalFinal - Xexp).^2)/(5*4-5))

% 殘差圖-----------------------------------------------------------
figure(1);
clf;
hold on;
plot(XcalFinal(:,1),Xexp(:,1),'m*');
plot(XcalFinal(:,2),Xexp(:,2),'r*');
plot(XcalFinal(:,3),Xexp(:,3),'b*');
plot(XcalFinal(:,4),Xexp(:,4),'k*');
legend('尾油','蠟油','輕油','氣+焦'); % 圖形注解
xlabel('Xcal,%'); % x軸注解
ylabel('Xexp,%'); % y軸注解
% title('Xexp-Xcal'); % 圖形標題
hold on;
x=0:0.1:80;
plot(x,x,'k-');

% y-t圖形顯示
figure(2);
clf;
hold on;
tf=max(time);
tspan=[0:1:tf];
[t,XcalFinal]=ode45(@modelEquation,tspan,X0,[],K,tf);
plot(t,XcalFinal(:,1),'m');
plot(t,XcalFinal(:,2),'r');
plot(t,XcalFinal(:,3),'b');
plot(t,XcalFinal(:,4),'k');
legend('尾油','蠟油','輕油','氣+焦'); % 圖形注解
xlabel('t/min'); % x軸注解
ylabel('X/%'); % y軸注解
% title('X--t'); % 圖形標題
hold;

%繪制試驗點
figure(2);
hold on;
plot(fourLumpData(:,1)',fourLumpData(:,2),'m*');
plot(fourLumpData(:,1)',fourLumpData(:,3),'r*');
plot(fourLumpData(:,1)',fourLumpData(:,4),'b*');
plot(fourLumpData(:,1)',fourLumpData(:,5),'k*');

% -----------------------------------------------------------------
function F = objFun(K,X0,Xexp,time)
tf=max(time);
tspan=[0:1:tf];
[nr,nc]=size(Xexp);
XcalFinal=zeros(nr,nc);
for i=1:nr
[t,Xcal]=ode45(@modelEquation,tspan,X0,[],K,tf);
n=find(t>=time(i));
m=n(1);
%線性內插
Xc=Xcal(m,

-(Xcal(m,

-Xcal(m-1,

)/(t(m)-t(m-1))*(t(m)-time(i));
XcalFinal(i,

=Xc;
end
objX = abs((XcalFinal - Xexp)./100).^2;
F=sum(sum(objX))^0.5;
error=(XcalFinal - Xexp)
avrError=sum(abs(error))/nr
avrErr=sum(avrError)/nc
relatError=(XcalFinal - Xexp)./Xexp*100
avrRelError=sum(abs(relatError))/nr
avrRelErr=sum(avrRelError)/nc
K

function dxdt = modelEquation(t,X,K,tf)
f1=-(K(1)+K(2)+K(3))*X(1);
f2= K(1)*X(1)-(K(4)+K(5))*X(2);
f3= K(2)*X(1)+K(4)*X(2);
f4= K(3)*X(1)+K(5)*X(2);
dxdt = [f1;f2;f3;f4];

贊一下

回復此樓

8樓2014-09-04 20:20:33

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

查看全部 10 個回答

月只藍

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

專家經驗: +1059
計算強帖: 8
應助: 1712 (講師)
貴賓: 8.888
金幣: 68127.7
散金: 1938
紅花: 443
沙發(fā): 4
帖子: 4373
在線: 3291.4小時
蟲號: 1122189
注冊: 2010-10-14
專業(yè): 宇宙學
管轄: 計算模擬區(qū)

【答案】應助回帖

感謝參與，應助指數 +1

參數估計問題或者叫常微分方程組參數擬合。參數太多，MATLAB求解估計困難不小，而且數據僅有4組。
建議你找有高版本1stOpt軟件的專家?guī)湍阕鲆幌拢扑]本版塊的專家dingd。

贊一下

回復此樓

MATLAB、MS小問題、普通問題請發(fā)帖求助！時間精力有限，恕不接受無償私信求助。

2樓2014-09-02 14:26:53

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

Cle_key

鐵蟲 (小有名氣)

應助: 22 (小學生)
金幣: 1560.8
紅花: 2
帖子: 178
在線: 28.9小時
蟲號: 2625972
注冊: 2013-08-31
性別: GG
專業(yè): 內流流體力學

【答案】應助回帖

感謝參與，應助指數 +1

雖然樓主自稱小女子。。。但肯定不幫

[ 發(fā)自手機版 http://www.gaoyang168.com/3g ]

贊一下

回復此樓

堅持就能成功

3樓2014-09-02 17:38:33

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

polypro

木蟲 (正式寫手)

應助: 197 (高中生)
金幣: 3921.3
散金: 160
紅花: 22
帖子: 512
在線: 62.6小時
蟲號: 1809438
注冊: 2012-05-11
性別: GG
專業(yè): 理論和計算化學

【答案】應助回帖

感謝參與，應助指數 +1

按自由約束條件，具有唯一解的情況是自由變量數=組分數-平衡相數+2，LZ的當前數據應該有很多不定解

贊一下

回復此樓

泉涸,魚相與處于陸,相呴以濕,相濡以沫,不如相忘于江湖。

4樓2014-09-03 08:20:26

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

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

[考研] 286求調劑 +6	lemonzzn 2026-03-16	9/450	2026-03-18 15:31 by vgtyfty
[考研] 0854可跨調劑，一作一項核心論文五項專利，省、國級證書40+數一英一287 +8	小李0854 2026-03-16	8/400	2026-03-18 14:35 by 搏擊518
[考研] 0817 化學工程 299分求調劑有科研經歷有二區(qū)文章 +7	rare12345 2026-03-18	7/350	2026-03-18 14:31 by laoshidan
[考研] 材料專碩306英一數二 +10	z1z2z3879 2026-03-16	13/650	2026-03-18 14:20 by 007_lilei
[考研] 265求調劑 +3	梁梁校校 2026-03-17	3/150	2026-03-18 09:12 by zhukairuo
[考研] 328求調劑，英語六級551，有科研經歷 +3	生物工程調劑 2026-03-16	8/400	2026-03-17 19:03 by Wangjingyue
[考研] 308求調劑 +4	是Lupa啊 2026-03-16	4/200	2026-03-17 17:12 by ruiyingmiao
[考研] 材料工程專碩274一志愿211求調劑 +6	薛云鵬 2026-03-15	6/300	2026-03-17 11:05 by 學員h26Tkc
[考研] 一志愿南京大學，080500材料科學與工程，調劑 +4	Jy? 2026-03-16	4/200	2026-03-17 11:02 by gaoqiong
[考研] 085600材料與化工求調劑 +13	enenenhui 2026-03-13	14/700	2026-03-16 15:19 by 了了了了。。
[考研] 285求調劑 +6	ytter 2026-03-12	6/300	2026-03-16 15:05 by njzyff
[考研] 中科大材料專碩319求調劑 +3	孟鑫材料 2026-03-13	3/150	2026-03-14 18:10 by houyaoxu
[考研] 297一志愿上交085600求調劑 +5	指尖八千里 2026-03-14	5/250	2026-03-14 17:26 by a不易
[考研] 297求調劑 +4	學海漂泊 2026-03-13	4/200	2026-03-14 11:51 by 熱情沙漠
[考研] 學碩285求調劑 +13	Wisjxn 2026-03-12	46/2300	2026-03-14 10:33 by JourneyLucky
[考研] 材料與化工求調劑一志愿 985 總分 295 +8	dream…… 2026-03-12	8/400	2026-03-13 22:17 by 星空星月
[考研] 281求調劑 +9	Koxui 2026-03-12	11/550	2026-03-13 20:50 by Koxui
[考研] 311求調劑 +3	冬十三 2026-03-13	3/150	2026-03-13 20:41 by JourneyLucky
[考研] 材料301分求調劑 +5	Liyouyumairs 2026-03-12	5/250	2026-03-13 14:42 by JourneyLucky
[考博] 26讀博 +4	Rui135246 2026-03-12	10/500	2026-03-13 07:15 by gaobiao