| 8 | 1/1 | 返回列表 |
| 查看: 2881 | 回復(fù): 7 | |||
ronghengyidu木蟲 (小有名氣)
|
[求助]
matlab擬合反應(yīng)動(dòng)力學(xué)參數(shù) 已有1人參與
|
|
本人是matlab大白,最近仿照黃華江編著的“實(shí)用化工計(jì)算機(jī)模擬”例題7-5擬合動(dòng)力學(xué)參數(shù), 動(dòng)力學(xué)方程為dc1/dt=k1*c1/(k2+k3*c1+k4*c2+k5*c1*c2),擬合k1,k2,k3,k4,k5, matlab程序和動(dòng)力學(xué)數(shù)據(jù)見附件,運(yùn)行出現(xiàn)如下錯(cuò)誤,請各位大拿幫助。 Error: File: KineticsEst.m Line: 47 Column: 60 Unbalanced or unexpected parenthesis or bracket. Error in run (line 64) evalin('caller', [script ';']); Error: File: KineticsEst.m Line: 47 Column: 60 Unbalanced or unexpected parenthesis or bracket. Error in run (line 64) evalin('caller', [script ';']); |
程序 |
|
樓主的問題,沒有看明白,找到了 黃華江編著的“實(shí)用化工計(jì)算機(jī)模擬”例題7-5擬合動(dòng)力學(xué)參數(shù) 7個(gè)微分方程,5個(gè)擬合參數(shù)k1,k2,k3,k4,k5,初值t=0時(shí),x為0.1883, 0.2507, 0.0467, 0.0899, 0.1804, 0.1394, 0.1046 q = 8.75 + k5, dx1/t= k5-q*x1- k1*x1*x2-k4*x1*x6*sqrt(0.9), dx2/t= 7.0-q*x2 - k1*x1*x2-2*k2*x2*x3, dx3/t= 1.75 -q*x3 - k2*x2*x3, dx4/t= -q*x4 + 2*k1*x1*x2-k3*x4*x5, dx5/t= -q*x5 + 3*k2*x2*x3-k3*x4*x5, dx6/t= -q*x6 + 2*k3*x4*x5-k4*x1*x6*sqrt(0.9), dx7/t= -q*x7 + 2*k4*x1*x6*sqrt(0.9) 數(shù)據(jù):t,x1,x4,x5,x6 0 0.1883 0.0899 0.1804 0.1394 0.0100 0.2047 0.0866 0.1729 0.1297 0.0200 0.2181 0.0856 0.1680 0.1205 0.0300 0.2291 0.0863 0.1647 0.1123 0.0400 0.2382 0.0878 0.1623 0.1053 0.0500 0.2459 0.0899 0.1604 0.0995 0.0600 0.2523 0.0921 0.1588 0.0948 0.0700 0.2576 0.0945 0.1574 0.0911 0.0800 0.2622 0.0968 0.1561 0.0882 0.0900 0.2660 0.0989 0.1548 0.0859 0.1000 0.2692 0.1010 0.1537 0.0842 0.1100 0.2719 0.1028 0.1525 0.0830 0.1200 0.2742 0.1045 0.1515 0.0821 0.1300 0.2761 0.1060 0.1505 0.0814 0.1400 0.2777 0.1074 0.1495 0.0810 0.1500 0.2790 0.1086 0.1487 0.0807 0.1600 0.2801 0.1096 0.1479 0.0805 0.1700 0.2811 0.1106 0.1471 0.0803 0.1800 0.2819 0.1114 0.1465 0.0803 0.1900 0.2825 0.1121 0.1458 0.0803 0.2000 0.2830 0.1127 0.1453 0.0803 用OpenLu求解: 結(jié)果(k1,k2,k3,k4,k5,目標(biāo)函數(shù)值): 17.60849027475542 73.06467445464085 51.32551162517818 23.02494410245102 6.001239504989786 6.794401934343191e-008 優(yōu)于原書matlab結(jié)果: 17.49, 72.2996, 50.9086, 22.5339, 5.9799, 3.614535245714516e-007 |
木蟲 (小有名氣)
送紅花一朵 |
非常感謝您的回復(fù),我只有一個(gè)微分方程,dc1/dt=k1*c1/(k2++k3*c1+k4*c2+k5*c1*c2), 擬合5個(gè)參數(shù)k1,k2,k3,k4,k5。我的t,c1,c2數(shù)據(jù)在附件中,matlab程序也在附件中,這兩個(gè)文件都可以用記事本打開。在這里我把程序和實(shí)驗(yàn)數(shù)據(jù)再貼出來。 function KineticsEst % 動(dòng)力學(xué)ODE方程模型的參數(shù)估計(jì) clear all clc k0 = [0.5 0.5 0.5 0.5 0.5]; % 參數(shù)初值 lb = [0 0 0 0 0]; % 參數(shù)下限 ub = [+inf +inf +inf +inf +inf]; % 參數(shù)上限 x0 = [150 1500]; KineticsData; yexp = ExpData(:,2:3); % yexp: 實(shí)驗(yàn)數(shù)據(jù)[x1 x2] % 使用函數(shù)fmincon()進(jìn)行參數(shù)估計(jì) [k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],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('\tk5 = %.4f\n',k(5)) 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 = ExpData(:,1); [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1:2) = x(:,1:2); f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2)); % ------------------------------------------------------------------ function f = ObjFunc4LNL(k,x0,yexp) tspan = = ExpData(:,1); [t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1:2) = x(:,1:2); f1 = y(:,1) - yexp(:,1); f2 = y(:,2) - yexp(:,2); f = [f1; f2]; % ------------------------------------------------------------------ function dxdt = KineticEqs(t,x,k) dxdt = ... k(1)*x(1)./(k(2)+k(3)*x(1)+k(4)*x(2)+k(5)*x(1).*x(2)); 實(shí)驗(yàn)數(shù)據(jù)如下: 動(dòng)力學(xué)數(shù)據(jù): t c1 c2 0 150 1500 10 137.3 1487.3 30 107.7 1457.7 60 72.8 1422.8 90 48.8 1398.8 120 22.9 1372.9 150 11.8 1361.8 180 5.4 1355.4 220 1.9 1351.9 240 0.9 1350.9 期望樓上能夠指點(diǎn),如果您能夠給出openLU程序,我也非常感謝,因?yàn)槲业谝淮温犝f這個(gè)程序,怕仿照編程又要出錯(cuò)。非常感謝您在節(jié)假日能夠回答我的問題。 |
|
是不是還有一個(gè)微分方程,dc2/dt=... ...。 否則應(yīng)該怎樣使用c2數(shù)據(jù)呢? |
木蟲 (小有名氣)
|
如果保留1500的c2那一列數(shù)據(jù),則微分方程為: dc1/dt=c1/(k1+k2*c1+k3*c1^2), dc2/dt=dc1/dt OpenLu代碼: 結(jié)果(k1,k2,k3,目標(biāo)函數(shù)值): -28.81145245946761 -0.6456992022469275 1.112333606324668e-003 39.20358725345449 可以看出,擬合參數(shù)值是相同的。 以下是c1、c2的計(jì)算值: 0. 150. 1500. 10. 135.485 1485.48 30. 108.705 1458.7 60. 73.9603 1423.96 90. 45.8941 1395.89 120. 25.0957 1375.1 150. 11.8172 1361.82 180. 4.86412 1354.86 220. 1.31349 1351.31 240. 0.665649 1350.67 |
| 8 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 268求調(diào)劑 +8 | 簡單點(diǎn)0 2026-03-17 | 8/400 |
|
|---|---|---|---|---|
|
[考研] 初試 317 +3 | 半拉月丙 2026-03-20 | 3/150 |
|
|
[考研] 能源材料化學(xué)課題組招收碩士研究生8-10名 +5 | 脫穎而出 2026-03-16 | 15/750 |
|
|
[考研] 一志愿重慶大學(xué)085700資源與環(huán)境總分308求調(diào)劑 +6 | 墨墨漠 2026-03-20 | 6/300 |
|
|
[考研] 301求調(diào)劑 +10 | yy要上岸呀 2026-03-17 | 10/500 |
|
|
[考研] 材料 336 求調(diào)劑 +3 | An@. 2026-03-18 | 4/200 |
|
|
[考研] 藥學(xué)383 求調(diào)劑 +3 | 藥學(xué)chy 2026-03-15 | 5/250 |
|
|
[考研] 求調(diào)劑一志愿南京航空航天大學(xué)289分 +3 | @taotao 2026-03-19 | 3/150 |
|
|
[考研] 一志愿西南交通 專碩 材料355 本科雙非 求調(diào)劑 +5 | 西南交通專材355 2026-03-19 | 5/250 |
|
|
[考研] 工科材料085601 279求調(diào)劑 +7 | 困于星晨 2026-03-17 | 9/450 |
|
|
[考研] 材料學(xué)碩318求調(diào)劑 +5 | February_Feb 2026-03-19 | 5/250 |
|
|
[考研] 一志愿中國海洋大學(xué),生物學(xué),301分,求調(diào)劑 +5 | 1孫悟空 2026-03-17 | 6/300 |
|
|
[考研] 0703化學(xué)調(diào)劑 +10 | 妮妮ninicgb 2026-03-15 | 14/700 |
|
|
[考研] 0703化學(xué)調(diào)劑 +5 | pupcoco 2026-03-17 | 8/400 |
|
|
[考研] 材料專碩306英一數(shù)二 +10 | z1z2z3879 2026-03-16 | 13/650 |
|
|
[考研] 0703化學(xué)調(diào)劑 +3 | 妮妮ninicgb 2026-03-17 | 3/150 |
|
|
[考研] 293求調(diào)劑 +11 | zjl的號 2026-03-16 | 16/800 |
|
|
[考研] 302求調(diào)劑 +4 | 小賈同學(xué)123 2026-03-15 | 8/400 |
|
|
[考研] 318求調(diào)劑 +3 | Yanyali 2026-03-15 | 3/150 |
|
|
[考研] 070305求調(diào)劑 +3 | mlpqaz03 2026-03-14 | 4/200 |
|