| 6 | 1/1 | 返回列表 |
| 查看: 1850 | 回復: 5 | ||
birdcanfly金蟲 (正式寫手)
|
[求助]
微分方程ode45求解,最小二乘法優(yōu)化微分方程參數(shù),程序運行求助 已有1人參與
|
|
仿照下面的帖子(http://www.gaoyang168.com/bbs/viewthread.php?tid=6425538&authorid=1122189)進行了程修改,求解微分方程組的系數(shù),微分方程組中包含了指前系數(shù)和活化能。 現(xiàn)在運行程序,matlab跑起來沒完,無結果,無報錯。請求高手指點,協(xié)助調通程序,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) % 目標函數(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))];%形成目標優(yōu)化函數(shù) function dCdt = KineticsEqs(t,C,k) % ODE模型方程,C為濃度,k為微分方程中的系數(shù),t為導數(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];%微分方程組 |
金蟲 (正式寫手)
|
可以順利運行,求出了參數(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) % 目標函數(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))];%形成目標優(yōu)化函數(shù) function dCdt = KineticsEqs(t,C,k) % ODE模型方程,C為濃度,k為微分方程中的系數(shù),t為導數(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];%微分方程組 |
金蟲 (正式寫手)
|
將參數(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 |
銅蟲 (正式寫手)
金蟲 (正式寫手)
| 6 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 317求調劑 +3 | 申子申申 2026-03-19 | 6/300 |
|
|---|---|---|---|---|
|
[考研] 286分人工智能專業(yè)請求調劑愿意跨考! +3 | lemonzzn 2026-03-17 | 3/150 |
|
|
[論文投稿]
申請回稿延期一個月,編輯同意了。但系統(tǒng)上的時間沒變,給編輯又寫郵件了,沒回復
10+3
|
wangf9518 2026-03-17 | 4/200 |
|
|
[考研] 一志愿蘇州大學材料求調劑,總分315(英一) +3 | sbdksD 2026-03-19 | 3/150 |
|
|
[考博] 申博26年 +3 | 八6八68 2026-03-19 | 3/150 |
|
|
[考研] 0703化學調劑 ,六級已過,有科研經歷 +12 | 曦熙兮 2026-03-15 | 12/600 |
|
|
[考研] 0703化學調劑 +5 | pupcoco 2026-03-17 | 8/400 |
|
|
[考研] 一志愿南昌大學,327分,材料與化工085600 +3 | Ncdx123456 2026-03-19 | 3/150 |
|
|
[考研] 一志愿中海洋材料工程專碩330分求調劑 +7 | 小材化本科 2026-03-18 | 7/350 |
|
|
[考研] 311求調劑 +6 | 26研0 2026-03-15 | 6/300 |
|
|
[考研] 298-一志愿中國農業(yè)大學-求調劑 +7 | 手機用戶 2026-03-17 | 7/350 |
|
|
[考研] 334求調劑 +3 | 志存高遠意在機?/a> 2026-03-16 | 3/150 |
|
|
[考研] 308求調劑 +4 | 是Lupa啊 2026-03-16 | 4/200 |
|
|
[考研] 一志愿蘇州大學材料工程(085601)專碩有科研經歷三項國獎兩個實用型專利一項省級立項 +6 | 大火山小火山 2026-03-16 | 8/400 |
|
|
[考博] 26申博 +4 | 八6八68 2026-03-16 | 4/200 |
|
|
[考研] 材料工程專碩274一志愿211求調劑 +6 | 薛云鵬 2026-03-15 | 6/300 |
|
|
[考研] 283求調劑 +3 | 聽風就是雨; 2026-03-16 | 3/150 |
|
|
[考研] 中科大材料與化工319求調劑 +3 | 孟鑫材料 2026-03-14 | 3/150 |
|
|
[考研] 復試調劑 +3 | 呼呼?~+123456 2026-03-14 | 3/150 |
|
|
[考研] 266求調劑 +4 | 學員97LZgn 2026-03-13 | 4/200 |
|