| 10 | 1/1 | 返回列表 |
| 查看: 882 | 回復(fù): 9 | ||
zhaoshazhu新蟲 (小有名氣)
|
[求助]
求Matlab指導(dǎo) 已有2人參與
|
| 為什么有時(shí)候Matlab給初值后有時(shí)計(jì)算快有時(shí)計(jì)算很慢,一直busy狀態(tài)?是出現(xiàn)死循環(huán)還是程序有問題?請求幫忙。 |
木蟲 (知名作家)
金蟲 (小有名氣)
新蟲 (小有名氣)
|
謝謝您的指導(dǎo),我現(xiàn)在就是初值選取有困難,很多時(shí)候都是自己嘗試,但是計(jì)算的結(jié)果不好,所得到參數(shù)的置信區(qū)間很大,請問您有好的建議嗎?下面是我做動(dòng)力學(xué)參數(shù)所編程的代碼。可以幫我看一下有沒有問題嗎? function DMMwan180 clear all clc k0 = [1 1 1 1 3 1 1 1 0.3 0.007]; lb = [0 0 0 0 0 0 0 0 0 0]; ub = [inf inf inf inf inf inf inf inf inf inf]; P0 =[0.015481 4.644213 0.340306 0 0 0; 0.015481 4.644213 0.340306 0 0 0; 0.018316 4.579047 0.402637 0 0 0; 0.022423 4.484655 0.492922 0 0 0; 0.018316 4.579047 0.402637 0 0 0; 0.022423 4.484655 0.492922 0 0 0; 0.040656 4.065616 0.893728 0 0 0; 0.015481 4.644213 0.340306 0 0 0; 0.015481 4.644213 0.340306 0 0 0; 0.018316 4.579047 0.402637 0 0 0; 0.022423 4.484655 0.492922 0 0 0; 0.040656 4.065616 0.893728 0 0 0 ]; % 初始分壓,MPa Pi=[0.0000439 4.6409577 0.3473380 0.0001568 0.0058211 0.0056825; 0.0000709 4.6414433 0.3484888 0.0003922 0.0067384 0.0028664; 0.0000841 4.5751124 0.4132686 0.0005222 0.0086355 0.0023772; 0.0002438 4.4758947 0.5083807 0.0014471 0.0118155 0.0022182; 0.0003697 4.5710387 0.4156405 0.0020687 0.0094840 0.0013984; 0.0009699 4.4759239 0.5083459 0.0032029 0.0100967 0.0014606; 0.0036297 4.0578876 0.9185219 0.0065201 0.0116013 0.0018393; 0.0001993 4.6368235 0.3512604 0.0010305 0.0088804 0.0018058; 0.0009229 4.6361031 0.3516653 0.0023232 0.0080420 0.0009435; 0.0015467 4.5720515 0.4148248 0.0039501 0.0068821 0.0007448; 0.0036621 4.4734942 0.5082864 0.0048581 0.0086486 0.0010507; 0.0087555 4.0513963 0.9187905 0.0087188 0.0105949 0.0017440 ]; % 經(jīng)過Wc/F0后,各物質(zhì)分壓,MPa % 使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì) opt=optimset('Algorithm','levenberg-marquardt'); [k,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@ObjFunc,k0,[],[],opt,P0,Pi); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\n') fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1)) fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2)) fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3)) fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4)) fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5)) fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6)) fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7)) fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8)) fprintf('\tk9 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9)) fprintf('\tk10 = %.4f ± %.4f\n',k(10),ci(10,2)-k(10)) fprintf('\t殘差平方和 = %.4f\n',resnorm) fprintf('\texitflag = %.4f\n',exitflag) fprintf('\tresidual = %.4f\n',residual) % ------------------------------------------------------------------ function f = ObjFunc(k,P0,Pi) % 目標(biāo)函數(shù) [m,n] = size(P0); Pcal = zeros(m,n); tspan = [0 2.05; 0 1.37; 0 1.62; 0 1.98; 0 1.21; 0 1.49; 0 2.70; 0 1.02; 0 0.82; 0 0.97; 0 1.19; 0 2.16 ]; % 即Wc/F0,g.h/mol for i = 1:m [t PP] = ode45(@Euqations,tspan(i, ,P0(i, ,[],k);Pcal(i, = PP(end, ;end f= Pcal-Pi; % ------------------------------------------------------------------ function dPdt = Euqations(t, P, k) % here t = Wc / F0 denom = 1+k(4)*P(1)+k(5)*P(3)+k(6)*P(4)+k(7)*P(5)+k(8)*P(6); % k(4) = KDMM,k(5) = KME ,k(6)=KHPM,k(7)=KPDO,k(8)=KNPA,k(9)=Kp1,k(10)=Kp2 theA =k(4)*(P(1)*P(2)-P(4)*P(3)/k(9)*P(2)) / denom^2; theB =k(6)*( P(4)*P(2)-P(5)*P(3)/k(10)*P(2))/ denom^2; theC =k(7)*P(5)*P(2)/denom^2; r1 = k(1)*theA; r2 = k(2)*theB; r3 = k(3)*theC; dPDMMdt = -r1; dPHdt = -2*r1-2*r2-r3; dPMEdt = r1+r2; dPHPMdt = r1-r2; dPPDOdt = r2-r3; dPNPAdt = r3; dPdt = [dPDMMdt;dPHdt;dPMEdt;dPHPMdt;dPPDOdt;dPNPAdt]; |
新蟲 (小有名氣)
木蟲 (知名作家)
新蟲 (小有名氣)
新蟲 (小有名氣)
|
這是我的算法您可以參考一下,希望得到您的幫助。 function DMMwan180 clear all clc k0 = [10 1 1 1 3 10 1 1 0.3 0.007]; lb = [0 0 0 0 0 0 0 0 0 0]; ub = [inf inf inf inf inf inf inf inf inf inf]; P0 =[0.015481 4.644213 0.340306 0 0 0; 0.015481 4.644213 0.340306 0 0 0; 0.018316 4.579047 0.402637 0 0 0; 0.022423 4.484655 0.492922 0 0 0; 0.018316 4.579047 0.402637 0 0 0; 0.022423 4.484655 0.492922 0 0 0; 0.040656 4.065616 0.893728 0 0 0; 0.015481 4.644213 0.340306 0 0 0; 0.015481 4.644213 0.340306 0 0 0; 0.018316 4.579047 0.402637 0 0 0; 0.022423 4.484655 0.492922 0 0 0; 0.040656 4.065616 0.893728 0 0 0 ]; % 初始分壓,MPa Pi=[0.0000439 4.6409577 0.3473380 0.0001568 0.0058211 0.0056825; 0.0000709 4.6414433 0.3484888 0.0003922 0.0067384 0.0028664; 0.0000841 4.5751124 0.4132686 0.0005222 0.0086355 0.0023772; 0.0002438 4.4758947 0.5083807 0.0014471 0.0118155 0.0022182; 0.0003697 4.5710387 0.4156405 0.0020687 0.0094840 0.0013984; 0.0009699 4.4759239 0.5083459 0.0032029 0.0100967 0.0014606; 0.0036297 4.0578876 0.9185219 0.0065201 0.0116013 0.0018393; 0.0001993 4.6368235 0.3512604 0.0010305 0.0088804 0.0018058; 0.0009229 4.6361031 0.3516653 0.0023232 0.0080420 0.0009435; 0.0015467 4.5720515 0.4148248 0.0039501 0.0068821 0.0007448; 0.0036621 4.4734942 0.5082864 0.0048581 0.0086486 0.0010507; 0.0087555 4.0513963 0.9187905 0.0087188 0.0105949 0.0017440 ]; % 經(jīng)過Wc/F0后,各物質(zhì)分壓,MPa % 使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì) opt=optimset('Algorithm','levenberg-marquardt'); [k,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@ObjFunc,k0,[],[],opt,P0,Pi); ci = nlparci(k,residual,jacobian); fprintf('\n\n使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\n') fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1)) fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2)) fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3)) fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4)) fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5)) fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6)) fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7)) fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8)) fprintf('\tk9 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9)) fprintf('\tk10 = %.4f ± %.4f\n',k(10),ci(10,2)-k(10)) fprintf('\t殘差平方和 = %.4f\n',resnorm) fprintf('\texitflag = %.4f\n',exitflag) fprintf('\tresidual = %.4f\n',residual) % ------------------------------------------------------------------ function f = ObjFunc(k,P0,Pi) % 目標(biāo)函數(shù) [m,n] = size(P0); Pcal = zeros(m,n); tspan = [0 2.05; 0 1.37; 0 1.62; 0 1.98; 0 1.21; 0 1.49; 0 2.70; 0 1.02; 0 0.82; 0 0.97; 0 1.19; 0 2.16 ]; % 即Wc/F0,g.h/mol for i = 1:m [t PP] = ode45(@Euqations,tspan(i, ,P0(i, ,[],k);Pcal(i, = PP(end, ;end f= Pcal-Pi; % ------------------------------------------------------------------ function dPdt = Euqations(t, P, k) % here t = Wc / F0 denom = 1+k(4)*P(1)+k(5)*P(3)+k(6)*P(4)+k(7)*P(5)+k(8)*P(6); % k(4) = KDMM,k(5) = KME ,k(6)=KHPM,k(7)=KPDO,k(8)=KNPA,k(9)=Kp1,k(10)=Kp2 theA =k(4)*(P(1)*P(2)-P(4)*P(3)/k(9)*P(2)) / denom^2; theB =k(6)*( P(4)*P(2)-P(5)*P(3)/k(10)*P(2))/ denom^2; theC =k(7)*P(5)*P(2)/denom^2; r1 = k(1)*theA; r2 = k(2)*theB; r3 = k(3)*theC; dPDMMdt = -r1; dPHdt = -2*r1-2*r2-r3; dPMEdt = r1+r2; dPHPMdt = r1-r2; dPPDOdt = r2-r3; dPNPAdt = r3; dPdt = [dPDMMdt;dPHdt;dPMEdt;dPHPMdt;dPPDOdt;dPNPAdt]; |
木蟲 (知名作家)
新蟲 (小有名氣)
| 10 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 求調(diào)劑 +5 | 爭取九點(diǎn)睡 2026-03-28 | 5/250 |
|
|---|---|---|---|---|
|
[考研] 材料與化工272求調(diào)劑 +7 | 阿斯蒂芬2004 2026-03-28 | 7/350 |
|
|
[考研] 一志愿南昌大學(xué)324求調(diào)劑 +7 | hanamiko 2026-03-27 | 7/350 |
|
|
[考研] 張芳銘-中國農(nóng)業(yè)大學(xué)-環(huán)境工程專碩-298 +4 | 手機(jī)用戶 2026-03-26 | 4/200 |
|
|
[考研] 086000調(diào)劑 +3 | 7901117076 2026-03-26 | 3/150 |
|
|
[考研] 化學(xué)調(diào)劑 +4 | 愛吃番茄的旭 2026-03-24 | 5/250 |
|
|
[考博] 26申博 +3 | 加油沖! 2026-03-26 | 3/150 |
|
|
[考研] 考研化學(xué)308分求調(diào)劑 +10 | 你好明天你好 2026-03-23 | 12/600 |
|
|
[考研] 298調(diào)劑 +3 | jiyingjie123 2026-03-27 | 3/150 |
|
|
[考研] 329求調(diào)劑 +7 | 鈕恩雪 2026-03-25 | 7/350 |
|
|
[考研] 349求調(diào)劑 +4 | 李木子啊哈哈 2026-03-25 | 4/200 |
|
|
[考研] 333求調(diào)劑 +7 | 87639 2026-03-21 | 12/600 |
|
|
[考研] 085600 材料與化工 329分求調(diào)劑 +9 | Mr. Z 2026-03-25 | 9/450 |
|
|
[考研] 299求調(diào)劑 +4 | 15188958825 2026-03-25 | 4/200 |
|
|
[考研] 290分調(diào)劑求助 +3 | 吉祥止止陳 2026-03-25 | 3/150 |
|
|
[考研] 300求調(diào)劑,材料科學(xué)英一數(shù)二 +5 | leaflight 2026-03-24 | 5/250 |
|
|
[考研] 277分求調(diào)劑,跨調(diào)材料 +3 | 考研調(diào)劑lxh 2026-03-24 | 3/150 |
|
|
[考研] 環(huán)境學(xué)碩288求調(diào)劑 +8 | 皮皮皮123456 2026-03-22 | 8/400 |
|
|
[論文投稿] 急發(fā)核心期刊論文 +3 | 賢達(dá)問津 2026-03-23 | 5/250 |
|
|
[考研] 293求調(diào)劑 +3 | 濤濤Wjt 2026-03-22 | 5/250 |
|