作者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