作者boffin (新的開始)
看板MATLAB
標題[問題] 解ODE中 希望負值變為0
時間Tue Apr 24 19:55:45 2012
有初始濃度與反應式
想要計算濃度隨時間的變化
但結果出來有幾筆為負值 可是濃度不能為負的
因為是用ode去解的
沒有頭緒該用什麼方法去做判斷
想請問大家有沒有什麼方式
先前用很小的dt下去跑 可是會跳出inf的問題
所以換成用ODE解
或是有沒有人做過類似的計算
感謝 :)
=== 附上code ==================================================
function ode1
vol=pi*2.6^2*15/1000; % L
UVtime=vol/2.5*60; % s
tspan = [0,UVtime];
y0 = [5.16e+18 0 0 0 7.889e+17 0 0 0 1.941e+19 2.457e+19];
options = odeset('RelTol',1e-15,'AbsTol',1e-15);
RH=99; % change here
T=298; % K
R=0.082; % atm.L.mol^-1.k^-1
N=6.02e23; % Avogadro constant
K1=2.6e-10; % s^-1
K2=6.0e-34*((T/300)^-2.3); % cm^6.mol^-1.s^-1
K3=4.1e-03; % s^-1
K4=8e-12*exp(-2060/T); % cm^3.mol^-1.s^-1
K5=2.2e-10; % cm^-3.mol^-1s^-1
K6=1.7e-12*exp(-940/T); % cm^3.mol^-1s^-1
K7=2.3e-13*exp(600/T)+4.3e-14*exp(1000/T); % cm^3.mol^-1
[t,y] = ode15s(@mm_rre,tspan,y0,options);
y1=y;
t1=t;
save y1data.mat y1 t1
function dy = mm_rre(t,y)
% MM_RRE Michaelis-Menten Reaction Rate Equation
dy = zeros(10,1);
dy(1) = -K1*y(1)+K3*y(8)+2*K4*(y(2)+y(4))*y(8)+K6*y(3)*y(8)+K7*y(6)*y(6);
dy(2) = 2*K1*y(1)-K2*y(2)*y(1)*y(10)-K4*(y(2)+y(4))*y(8);
dy(3) = K5*y(4)*y(5)-K6*y(3)*y(5);
dy(4) = K3*y(8)-K5*y(4)*y(5);
dy(5) = -K5*y(4)*y(5);
dy(6) = K6*y(3)*y(5)-K7*y(6)*y(6);
dy(7) = K7*y(6)*y(6);
dy(8) = K2*y(2)*y(1)*y(10)-K4*(y(2)+y(4))*y(8)-K6*y(3)*y(5);
dy(9) = 0;
dy(10)= -K2*y(2)*y(1)*y(10);
end
end
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.112.66.253
1F：推 parker16:改變數 x=e^(x_1)就不會有負的 04/25 13:58
2F：→ boffin:odeset中加入 'nonnegtive',[哪幾項] 就OK了! 05/08 00:26