 |
 |
 |
 |
 |
 |
||
 |
 |
 |
 |
 |
| 6 | 1/1 | 返回列表 |
| 查看: 1854 | 回復(fù): 5 | ||
birdcanfly金蟲 (正式寫手)
|
[求助]
微分方程ode45求解，最小二乘法優(yōu)化微分方程參數(shù)，程序運行求助 已有1人參與
|
|
仿照下面的帖子（http://www.gaoyang168.com/bbs/viewthread.php?tid=6425538&authorid=1122189）進行了程修改，求解微分方程組的系數(shù)，微分方程組中包含了指前系數(shù)和活化能。 現(xiàn)在運行程序，matlab跑起來沒完，無結(jié)果，無報錯。請求高手指點，協(xié)助調(diào)通程序，50金幣相謝！ *********************** 程序代碼如下： function k1k2k3 format long clear all clc k0 = [1e40 1e20 1e40 1e30 100e3 50e3 100e3 100e3 0.8]; lb = [0 0 0 0 0 0 0 0 0]; ub = [inf inf inf inf inf inf inf inf 1]; data=... [ 0.00 100.000 0.00000 0.00000 0.00000 0.21 78.8686 18.802 1.99667 0 0.22 59.2346 37.9368 2.66223 0.332779 0.23 48.0865 46.7554 4.49251 0.166389 0.24 37.1048 57.0715 5.32446 0.665557 0.25 29.6173 61.2313 8.31947 0.831947 0.253 31.1148 59.5674 9.65058 0.831947 0.258 29.7837 61.2313 7.8203 0.998336 0.26 22.629 63.3943 12.3128 1.16473 0.263 24.792 60.3993 13.4775 1.33111 0.266 21.9634 62.396 13.9767 1.4975 0.27 21.2978 63.228 13.1448 1.83028 0.272 24.6256 59.0682 14.3095 1.66389 0.277 17.3045 58.4027 21.1314 3.1614 0.28 17.4709 58.7354 20.9651 2.82862 0.29 23.7937 49.9168 22.4626 3.32779 0.3 18.3028 43.7604 31.1148 6.48918 ]; x0=data(1,2:end); tspan=data(:,1)'; yexp = [data(2:end,2) data(2:end,3) data(2:end,4) data(2:end,5)]; [k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函數(shù)lsqnonlin()估計得到的參數(shù)值為:\n') fprintf('\tk1 = %.9f ± %.9f\n',k(1),ci(1,2)-k(1)) fprintf('\tk2 = %.9f ± %.9f\n',k(2),ci(2,2)-k(2)) fprintf('\tk3 = %.9f ± %.9f\n',k(3),ci(3,2)-k(3)) fprintf('\tk4 = %.9f ± %.9f\n',k(4),ci(4,2)-k(4)) fprintf('\tEa1 = %.9f ± %.9f\n',k(5),ci(5,2)-k(5)) fprintf('\tEa2 = %.9f ± %.9f\n',k(6),ci(6,2)-k(6)) fprintf('\tEa3 = %.9f ± %.9f\n',k(7),ci(7,2)-k(7)) fprintf('\tEa4 = %.9f ± %.9f\n',k(8),ci(8,2)-k(8)) fprintf('\ta = %.9f ± %.9f\n',k(9),ci(9,2)-k(9)) fprintf(' The sum of the squares is: %.9e\n\n',resnorm) ts=0  (max(tspan)-min(tspan))/100):max(tspan);[ts ys] = ode45(@KineticsEqs,ts,x0,[],k); yy = [data(:,2) data(:,3) data(:,4)]; I100=100*ones(16,1); plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo'); hold on plot(ts,ys(:,2)+ys(:,3),'r',tspan,yy(:,2),'r*'); plot(ts,ys(:,4),'k',tspan,yy(:,3),'k+'); plot(ts,I100(:,1)-ys(:,1)-ys(:,2)-ys(:,3)-ys(:,4),'g',tspan,yy(:,4),'g<') legend('C1的計算值','C1的實驗值','C2的計算值*5','C2的實驗值*5','C3的計算值*5','C3的實驗值*5','C4的計算值','C4的實驗值') function f = ObjFunc(k,tspan,x0,yexp) % 目標(biāo)函數(shù) [t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);% 求解常微分方程，其中tspan為t的取值點，x0為微分方程組的初始值，k為微分方程的系數(shù)，返回t為微分方程組解的取值點，Xsim為微分方程組的解 Xsim1=Xsim(:,1);%提取Xsim的第一列 Xsim2=Xsim(:,2);%提取Xsim的第二列 Xsim3=Xsim(:,3);%提取Xsim的第三列 Xsim4=Xsim(:,4);%提取Xsim的第四列 ysim(:,1) = Xsim1(2:end);%微分方程的第一個變量，從第二個解到最后一個解賦值給ysim的第一列 ysim(:,2) = Xsim2(2:end);%微分方程的第二個變量，從第二個解到最后一個解賦值給ysim的第二列 ysim(:,3) = Xsim3(2:end);%微分方程的第三個變量，從第二個解到最后一個解賦值給ysim的第三列 ysim(:,4) = Xsim4(2:end);%微分方程的第四個變量，從第二個解到最后一個解賦值給ysim的第四列 I100=100*ones(16,1); f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)+ysim(:,3)-yexp(:,2)) (ysim(:,4)-yexp(:,3)) (I100-ysim(:,1)-ysim(:,2)-ysim(:,3)-ysim(:,4)-yexp(:,4))];%形成目標(biāo)優(yōu)化函數(shù) function dCdt = KineticsEqs(t,C,k) % ODE模型方程，C為濃度，k為微分方程中的系數(shù)，t為導(dǎo)數(shù) a=k(9); Ea1=k(5); Ea2=k(6); Ea3=k(7); Ea4=k(8); R=8.314; Ttime=(-960.5*t^2+842.1*t+56.66+273.15); dC1dt =-k(1)*exp(-Ea1/(R*Ttime))*C(1);%第一個微分方程，等號前面是C(1)對t的微分 dC2dt =k(1)*exp(-Ea1/(R*Ttime))*a*C(1)-k(2)*exp(-Ea2/(R*Ttime))*C(2);%第二個微分方程，等號前面是C(2)對t的微分 dC3dt =k(2)*exp(-Ea2/(R*Ttime))*C(2)-k(3)*exp(-Ea3/(R*Ttime))*C(3);%第三個微分方程，等號前面是C(3)對t的微分 dC4dt =k(3)*exp(-Ea3/(R*Ttime))*C(3)-k(4)*exp(-Ea4/(R*Ttime))*C(4);%第四個微分方程，等號前面是C(4)對t的微分 dCdt = [dC1dt; dC2dt;dC3dt;dC4dt];%微分方程組 |
回復(fù)此樓» 猜你喜歡
哈爾濱理工大學(xué)物理系招收物理學(xué)考研調(diào)劑
已經(jīng)有6人回復(fù)
-
0702一志愿吉大B區(qū)求調(diào)劑
已經(jīng)有5人回復(fù)
-
物理學(xué)I論文潤色/翻譯怎么收費?
已經(jīng)有126人回復(fù)
-
求調(diào)劑
已經(jīng)有0人回復(fù)
-
0702一志愿吉大B區(qū)求調(diào)劑有論文
已經(jīng)有0人回復(fù)
-
請問還有沒有用Latex寫文章的小伙伴們？
已經(jīng)有0人回復(fù)
-
光學(xué)工程學(xué)碩調(diào)劑信息
已經(jīng)有31人回復(fù)
-
歡迎加入課題組
已經(jīng)有0人回復(fù)
-
散金幣，求好運，祝面上順利！
已經(jīng)有65人回復(fù)
-
華東師范大學(xué)芯片設(shè)計徐瓏真實水平咋樣
已經(jīng)有1人回復(fù)
-
調(diào)劑
已經(jīng)有0人回復(fù)
» 本主題相關(guān)價值貼推薦，對您同樣有幫助:
- 有方程，怎么求多元回歸方程的參數(shù)
已經(jīng)有12人回復(fù)
- 求大神幫忙matlab擬合函數(shù)求參數(shù)
已經(jīng)有9人回復(fù)
- 如何求解帶有固定參數(shù)的超定方程組？
已經(jīng)有21人回復(fù)
- 微分方程組求參數(shù)問題，求高人指教，感謝
已經(jīng)有11人回復(fù)
- matlab求解方程中的參數(shù)
已經(jīng)有21人回復(fù)
- 動力學(xué)參數(shù)擬合
已經(jīng)有26人回復(fù)
- 一階微分方程求解，并優(yōu)化參數(shù)值
已經(jīng)有14人回復(fù)
- matlab擬合方程參數(shù)時初值的選擇
已經(jīng)有15人回復(fù)
- matlab解微分方程
已經(jīng)有10人回復(fù)
- 用matlab求解方程出問題，請幫忙看看
已經(jīng)有3人回復(fù)
- 求助利用matlab的ODE45求解微分方程
已經(jīng)有16人回復(fù)
- 用Matlab求解方程
已經(jīng)有6人回復(fù)
- 最小二乘法解方程組系數(shù)
已經(jīng)有10人回復(fù)
- 用MATLAB ode45求解2階微分方程
已經(jīng)有4人回復(fù)
- 【求助】用matlab最優(yōu)化方法進行參數(shù)擬合
已經(jīng)有17人回復(fù)
- 【求助】急求各位大俠，關(guān)于用ode45解微分方程
已經(jīng)有8人回復(fù)
- 【求助】如何計算實驗數(shù)據(jù)與已知方程的相關(guān)系數(shù)？
已經(jīng)有6人回復(fù)
- 【求助】最小二乘法編寫指數(shù)擬合方程
已經(jīng)有19人回復(fù)
- 【求助】使用Matlab預(yù)估動力學(xué)方程問題
已經(jīng)有13人回復(fù)
2樓2014-09-18 22:12:01
birdcanfly
金蟲 (正式寫手)
- 應(yīng)助: 8 (幼兒園)
- 金幣: 735.7
- 散金: 406
- 紅花: 3
- 帖子: 320
- 在線: 105.9小時
- 蟲號: 399057
- 注冊: 2007-06-11
- 性別: GG
- 專業(yè): 生物醫(yī)用高分子材料
|
可以順利運行，求出了參數(shù)，可以很好的擬合文獻中的數(shù)據(jù)，但是參數(shù)和文獻報道的參數(shù)差好幾個數(shù)量級。 ********************* function k1k2k3 format short e clear all clc k0 = [1.0 1.0 1.0 1.0 139 59 169 100 0.8]; lb = [0 0 0 0 0 0 0 0 0]; ub = [inf inf inf inf inf inf inf inf 1]; data=... [ 0.21 78.8686 18.802 1.99667 0 0.22 59.2346 37.9368 2.66223 0.332779 0.23 48.0865 46.7554 4.49251 0.166389 0.24 37.1048 57.0715 5.32446 0.665557 0.25 29.6173 61.2313 8.31947 0.831947 0.253 31.1148 59.5674 9.65058 0.831947 0.258 29.7837 61.2313 7.8203 0.998336 0.26 22.629 63.3943 12.3128 1.16473 0.263 24.792 60.3993 13.4775 1.33111 0.266 21.9634 62.396 13.9767 1.4975 0.27 21.2978 63.228 13.1448 1.83028 0.272 24.6256 59.0682 14.3095 1.66389 0.277 17.3045 58.4027 21.1314 3.1614 0.28 17.4709 58.7354 20.9651 2.82862 0.29 23.7937 49.9168 22.4626 3.32779 0.3 18.3028 43.7604 31.1148 6.48918 ]; x0=data(1,2:end); tspan=data(:,1); yexp = [data(2:end,2) data(2:end,3) data(2:end,4) data(2:end,5)]; [k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函數(shù)lsqnonlin()估計得到的參數(shù)值為:\n') fprintf('\tk1 = %e ± %e\n',k(1),ci(1,2)-k(1)) fprintf('\tk2 = %e ± %e\n',k(2),ci(2,2)-k(2)) fprintf('\tk3 = %e ± %e\n',k(3),ci(3,2)-k(3)) fprintf('\tk4 = %e ± %e\n',k(4),ci(4,2)-k(4)) fprintf('\tEa1 = %e ± %e\n',k(5),ci(5,2)-k(5)) fprintf('\tEa2 = %e ± %e\n',k(6),ci(6,2)-k(6)) fprintf('\tEa3 = %e ± %e\n',k(7),ci(7,2)-k(7)) fprintf('\tEa4 = %e ± %e\n',k(8),ci(8,2)-k(8)) fprintf('\ta = %f ± %f\n',k(9),ci(9,2)-k(9)) ts=0.210 (max(tspan)-min(tspan))/100):max(tspan);[ts ys] = ode45(@KineticsEqs,ts,x0,[],k); yy = [data(:,2) data(:,3) data(:,4) data(:,5)]; I100=100*ones(101,1); plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo'); hold on plot(ts,ys(:,2)+ys(:,3),'r',tspan,yy(:,2),'r*'); plot(ts,ys(:,4),'k',tspan,yy(:,3),'k+'); plot(ts,I100-ys(:,1)-ys(:,2)-ys(:,3)-ys(:,4),'g',tspan,yy(:,4),'g<') legend('C1的計算值','C1的實驗值','C2的計算值','C2的實驗值','C3的計算值','C3的實驗值','C4的計算值','C4的實驗值') function f = ObjFunc(k0,tspan,x0,yexp) % 目標(biāo)函數(shù) [t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k0);% 求解常微分方程，其中tspan為t的取值點，x0為微分方程組的初始值，k為微分方程的系數(shù)，返回t為微分方程組解的取值點，Xsim為微分方程組的解 Xsim1=Xsim(:,1);%提取Xsim的第一列 Xsim2=Xsim(:,2);%提取Xsim的第二列 Xsim3=Xsim(:,3);%提取Xsim的第三列 Xsim4=Xsim(:,4);%提取Xsim的第四列 ysim(:,1) = Xsim1(2:end);%微分方程的第一個變量，從第二個解到最后一個解賦值給ysim的第一列 ysim(:,2) = Xsim2(2:end);%微分方程的第二個變量，從第二個解到最后一個解賦值給ysim的第二列 ysim(:,3) = Xsim3(2:end);%微分方程的第三個變量，從第二個解到最后一個解賦值給ysim的第三列 ysim(:,4) = Xsim4(2:end);%微分方程的第四個變量，從第二個解到最后一個解賦值給ysim的第四列 I100=100*ones(15,1); f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)+ysim(:,3)-yexp(:,2)) (ysim(:,4)-yexp(:,3)) (I100-ysim(:,1)-ysim(:,2)-ysim(:,3)-ysim(:,4)-yexp(:,4))];%形成目標(biāo)優(yōu)化函數(shù) function dCdt = KineticsEqs(t,C,k) % ODE模型方程，C為濃度，k為微分方程中的系數(shù)，t為導(dǎo)數(shù) a=k(9); Ea1=k(5); Ea2=k(6); Ea3=k(7); Ea4=k(8); R=8.314; Ttime=(-960.5*t^2+842.1*t+56.66+273.15); dC1dt =-k(1)*exp(-Ea1/(R*Ttime))*C(1);%第一個微分方程，等號前面是C(1)對t的微分 dC2dt =k(1)*exp(-Ea1/(R*Ttime))*a*C(1)-k(2)*exp(-Ea2/(R*Ttime))*C(2);%第二個微分方程，等號前面是C(2)對t的微分 dC3dt =k(2)*exp(-Ea2/(R*Ttime))*C(2)-k(3)*exp(-Ea3/(R*Ttime))*C(3);%第三個微分方程，等號前面是C(3)對t的微分 dC4dt =k(3)*exp(-Ea3/(R*Ttime))*C(3)-k(4)*exp(-Ea4/(R*Ttime))*C(4);%第四個微分方程，等號前面是C(4)對t的微分 dCdt = [dC1dt; dC2dt;dC3dt;dC4dt];%微分方程組 |
3樓2014-09-19 13:17:46
birdcanfly
金蟲 (正式寫手)
- 應(yīng)助: 8 (幼兒園)
- 金幣: 735.7
- 散金: 406
- 紅花: 3
- 帖子: 320
- 在線: 105.9小時
- 蟲號: 399057
- 注冊: 2007-06-11
- 性別: GG
- 專業(yè): 生物醫(yī)用高分子材料
|
將參數(shù)的初始值設(shè)定為文獻中報道的數(shù)值，報錯如下： Local minimum possible. lsqnonlin stopped because the size of the current step is less than the default value of the step size tolerance. <stopping criteria details> Warning: Matrix is singular to working precision. > In nlparci at 104 In ODEParmafit_non_test at 33 |
4樓2014-09-19 13:22:20
whqs8426212
銅蟲 (正式寫手)
- 應(yīng)助: 5 (幼兒園)
- 金幣: 13.2
- 散金: 375
- 紅花: 7
- 帖子: 782
- 在線: 685.7小時
- 蟲號: 2572504
- 注冊: 2013-07-30
- 專業(yè): 數(shù)論
5樓2014-09-20 18:18:39
birdcanfly
金蟲 (正式寫手)
- 應(yīng)助: 8 (幼兒園)
- 金幣: 735.7
- 散金: 406
- 紅花: 3
- 帖子: 320
- 在線: 105.9小時
- 蟲號: 399057
- 注冊: 2007-06-11
- 性別: GG
- 專業(yè): 生物醫(yī)用高分子材料
6樓2014-09-23 18:08:55
相關(guān)版塊跳轉(zhuǎn)
我要訂閱樓主 birdcanfly 的主題更新
| 6 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 08工學(xué)調(diào)劑 +5 | 用戶573181 2026-03-20 | 5/250 |
|
|---|---|---|---|---|
|
[考研] 274求調(diào)劑 +8 | S.H1 2026-03-18 | 8/400 |
|
|
|
[考研] 工科材料085601 279求調(diào)劑 +6 | 困于星晨 2026-03-17 | 8/400 |
|
|
|
[考研] 286分人工智能專業(yè)請求調(diào)劑愿意跨考！ +3 | lemonzzn 2026-03-17 | 4/200 |
|
|
|
[考研] 【考研調(diào)劑】化學(xué)專業(yè) 281分，一志愿四川大學(xué)，誠心求調(diào)劑 +6 | 吃吃吃才有意義 2026-03-19 | 6/300 |
|
|
|
[考研] 286求調(diào)劑 +6 | lemonzzn 2026-03-16 | 10/500 |
|
|
|
[考研] 一志愿中海洋材料工程專碩330分求調(diào)劑 +7 | 小材化本科 2026-03-18 | 7/350 |
|
|
|
[考研] 一志愿華中科技大學(xué)，080502，354分求調(diào)劑 +4 | 守候夕陽CF 2026-03-18 | 4/200 |
|
|
|
[考研] 328求調(diào)劑，英語六級551，有科研經(jīng)歷 +3 | 生物工程調(diào)劑 2026-03-17 | 7/350 |
|
|
|
[考研] 一志愿武理材料305分求調(diào)劑 +5 | 想上岸的鯉魚 2026-03-18 | 6/300 |
|
|
|
[考研] 085601專碩，總分342求調(diào)劑，地區(qū)不限 +5 | share_joy 2026-03-16 | 5/250 |
|
|
|
[考研] 收復(fù)試調(diào)劑生 +4 | 雨后秋荷 2026-03-18 | 4/200 |
|
|
|
[考研] 303求調(diào)劑 +4 | 睿08 2026-03-17 | 6/300 |
|
|
|
[考研] 278求調(diào)劑 +5 | 煙火先于春 2026-03-17 | 5/250 |
|
|
|
[考研] 290求調(diào)劑 +3 | p asserby. 2026-03-15 | 4/200 |
|
|
|
[考研] 一志愿蘇州大學(xué)材料工程（085601）專碩有科研經(jīng)歷三項國獎兩個實用型專利一項省級立項 +6 | 大火山小火山 2026-03-16 | 8/400 |
|
|
|
[考博] 26申博 +4 | 八6八68 2026-03-16 | 4/200 |
|
|
|
[考研] [導(dǎo)師推薦]西南科技大學(xué)國防/材料導(dǎo)師推薦 +3 | 尖角小荷 2026-03-16 | 6/300 |
|
|
|
[考研] 中科院材料273求調(diào)劑 +4 | yzydy 2026-03-15 | 4/200 |
|
|
|
[考研] 求老師收留調(diào)劑 +4 | jiang姜66 2026-03-14 | 5/250 |
|