| 5 | 2/1 | 返回列表 |
| 查看: 3613 | 回復(fù): 15 | |||
| 當(dāng)前只顯示滿足指定條件的回帖,點擊這里查看本話題的所有回帖 | |||
[求助]
matlab擬合方程參數(shù)時初值的選擇
|
|||
|
A+B=C A+W=R R+B=C+W 先考慮的三個反應(yīng)都是可逆反應(yīng),有以下幾個式子:編寫程序如下 %30℃ clear all clc global keq1 keq2 t0 c0 keq1 = 3.94; % 平衡常數(shù)k1 keq3 = 1.08; % 平衡常數(shù)k3 t0 = [5,10, 15, 30, 45, 60,90,120,160,200,240,300,360]; c0 = [5.939601,5.203973,0,0.022484,0.078956] % A B C R W 組成的初始濃度 ci=[3.742508811 4.213566482 0.990406938 0.252498405 0.090415512; 3.045716838 3.387013322 1.816960097 0.215546741 0.081072022; 2.504754523 2.901991595 2.301981824 0.193427102 0.099462604; 1.65098702 2.036686653 3.167286766 0.170236582 0.095144044; 1.271633833 1.655460403 3.548513016 0.124730781 0.10334349; 1.129887499 1.630815394 3.573158025 0.129798695 0.097764543; 1.01239208 1.460748979 3.74322444 0.144649391 0.073614958; 1.007417603 1.429180292 3.774793127 0.123477124 0.094249307; 0.990005723 1.499216939 3.704756481 0.135524067 0.076072022; 0.970073726 1.493966669 3.71000675 0.130791851 0.076387812; 0.96903024 1.429705429 3.77426799 0.119510017 0.101202216; 0.980873385 1.457988124 3.745985296 0.124189274 0.097022161; 0.98477262 1.453402965 3.750570455 0.1200936 0.0762608]; % cA cC cR 對應(yīng)時間數(shù)值 k0=[1 1 1]; % k1+, k2, k3+ 初始向值 lb = [0 0 0]; ub = [1000 1000 500]; % 上下限 % 使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計 [k,resnorm,residual,exitflag] = lsqnonlin(@ObjFunc,k0,lb,ub,[],ci); k1plus=k(1); k1minus= k1plus/keq1; k2=k(2); k3plus=k(3); k3minus= k3plus/keq3; % ------------------------------------------------------------------ function f = ObjFunc(k,ci) % 目標(biāo)函數(shù) global t0 c0 [t,c_cal] = ode45(@Euqations,t0,c0,[],k); f= sum((c_cal-ci)^2); % ------------------------------------------------------------------ function dcdt = Euqations(t,c,k) global keq1 keq2 k1plus=k(1); k1minus= k1plus/keq1; k2=k(2); k3plus=k(3); k3minus= k3plus/keq3; cA=c(:,1); cB=c(:,2); cC=c(:,3); cR=c(:,4); cW=c(:,5); dcAdt = k1plus*cA*cB-k1minus*cC+k2*cA*cW; dcCdt = -k1plus*cA*cB+k1minus*cC+k3plus*cB*cR-k3minus*cC*cW; dcRdt = k2*cA*cW-k3plus*cB*cR+k3minus*cC*cW; dcdt = [dcAdt; dcCdt;dcRdt]; 但是k的初值一直影響結(jié)果,程序要怎么改一下初值對結(jié)果影響不大(這個程序有問題,) |
木蟲 (正式寫手)
|
function KineticsEst clear all clc x0 = [4.639975 5.289842 0 0.039974 0.063806]; % A B C R W 組成的初始濃 k0=[1 1 1]; % k1, k2, k3 初始向值 lb = [0 0 0]; ub = [1 +inf +inf]; % 上下限 KineticsData3; yexp = Kinetics(:,2:6); % 使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計 [k,resnorm,residual,exitflag,output,lambda,jacobian]=... lsqnonlin(@ObjFunc,k0,lb,ub,[],x0,yexp); fprintf('\n\n使用函數(shù)lsqnonlin()估計得到的參數(shù)值為:\n') fprintf('\tk1 = %.4f ± %.4f\n',k(1)) fprintf('\tk2 = %.4f ± %.4f\n',k(2)) fprintf('\tk3 = %.4f ± %.4f\n',k(3)) % ------------------------------------------------------------------ function f = ObjFunc(k,x0,yexp) % 目標(biāo)函數(shù) tspan = [0;5;10;15;30;46;60;90;120;160;200;240;300;360]; [t x] = ode45(@Euqations,tspan,x0,[],k); y(:,1) = x(:,1); y(:,2:5) = x(:,2:5); f1 = y(:,1) - yexp(:,1); f2 = y(:,2) - yexp(:,2); f3 = y(:,3) - yexp(:,3); f4 = y(:,4) - yexp(:,4); f5 = y(:,5) - yexp(:,5); f = [f1; f2; f3; f4;f5]; % ------------------------------------------------------------------ function dcdt = Euqations(t,x,k) r1 = k(1)*x(1)*x(2)-k(1)*x(3)/3.94; r2 = k(2)*x(1)*x(5); r3 = k(3)*x(2)*x(4)-k(3)*x(3)*x(5)/1.08; dcAdt = -r1-r2; dcBdt = -r1-r3; dcCdt = r1+r3; dcRdt = r2-r3; dcWdt = -r2+r3; dcdt = [dcAdt;dcBdt; dcCdt;dcRdt;dcWdt;]; |
鐵桿木蟲 (職業(yè)作家)
|
我寫在附件里面了 |
鐵桿木蟲 (職業(yè)作家)
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 280分求調(diào)劑 一志愿085802 +4 | PUMPT 2026-03-22 | 7/350 |
|
|---|---|---|---|---|
|
[考研] 石河子大學(xué)(211、雙一流)碩博研究生長期招生公告 +3 | 李子目 2026-03-22 | 3/150 |
|
|
[考研] 生物學(xué)071000 329分求調(diào)劑 +5 | 我愛生物生物愛?/a> 2026-03-17 | 5/250 |
|
|
[考研] 289求調(diào)劑 +7 | 懷瑾握瑜l 2026-03-20 | 7/350 |
|
|
[考研] 319求調(diào)劑 +4 | 小力氣珂珂 2026-03-20 | 4/200 |
|
|
[考研] 求調(diào)劑 +5 | Zhangbod 2026-03-21 | 7/350 |
|
|
[考博] 招收博士1-2人 +3 | QGZDSYS 2026-03-18 | 4/200 |
|
|
[考研]
|
Grand777 2026-03-21 | 3/150 |
|
|
[考研] 一志愿深大,0703化學(xué),總分302,求調(diào)劑 +4 | 七月-七七 2026-03-21 | 4/200 |
|
|
[考研] 278求調(diào)劑 +9 | 煙火先于春 2026-03-17 | 9/450 |
|
|
[考研] 0805材料320求調(diào)劑 +3 | 深海物語 2026-03-20 | 3/150 |
|
|
[考研] 材料學(xué)學(xué)碩080502 337求調(diào)劑-一志愿華中科技大學(xué) +4 | 順順順mr 2026-03-18 | 5/250 |
|
|
[考研] 一志愿西南交大,求調(diào)劑 +5 | 材化逐夢人 2026-03-18 | 5/250 |
|
|
[考研] 321求調(diào)劑 +9 | 何潤采123 2026-03-18 | 11/550 |
|
|
[考研] 求調(diào)劑 +3 | eation27 2026-03-20 | 3/150 |
|
|
[考研] 工科材料085601 279求調(diào)劑 +7 | 困于星晨 2026-03-17 | 9/450 |
|
|
[考研] 0703化學(xué)調(diào)劑 +4 | 18889395102 2026-03-18 | 4/200 |
|
|
[考研] 材料專碩306英一數(shù)二 +10 | z1z2z3879 2026-03-16 | 13/650 |
|
|
[碩博家園] 湖北工業(yè)大學(xué) 生命科學(xué)與健康學(xué)院-課題組招收2026級食品/生物方向碩士 +3 | 1喜春8 2026-03-17 | 5/250 |
|
|
[考研] [導(dǎo)師推薦]西南科技大學(xué)國防/材料導(dǎo)師推薦 +3 | 尖角小荷 2026-03-16 | 6/300 |
|