ODE45: doubts about the result. Correct or not?

4 次查看（过去 30 天）
Marco Sammito
Marco Sammito 2016-11-28
评论： bio lim 2016-11-28
Hi. I have to solve
I do not know if it may be helpful, but note that a1=2.488125⋅10^14 is very big while the timespan is very small (tf=70⋅10^−6). I get the result without MATLAB returning any error or warning, but this result does not convince me (I expected the peak to be around t=4⋅10^−5). I have already tried switching to ode15s, but the output is the same. What else can I do to see if the output changes or not and then be sure that this result is correct?
function vdot = no_thermal_effect_98(t,v)
rho = 959.78;
P0 = 101325;
Pvap = 94285.313;
a1 = (P0 - 1800) / (20e-6)^2; %nondimensional
a2 = 2 * (1800 - P0) / (20e-6); %nondimensional
%
vdot = zeros(2,1);
vdot(1) = v(2);
vdot(2) = -1.5 * v(2) * v(2) / v(1) + 1 / (v(1) * rho) *...
(Pvap - P0 - a1 * t^2 - a2 * t);
run file
x0 = 5e-6; %meters
tf = 70e-6; %seconds
[t,v] = ode45(@no_thermal_effect_98,[0,tf],[x0,0]);
[t,v(:,1)];
plot(t,v(:,1))
  1 个评论
bio lim
bio lim 2016-11-28
I think the main problem arises when
vdot(2) = -1.5 * v(2) * v(2) / v(1) + 1 / (v(1) * rho) *...
(Pvap - P0 - a1 * t^2 - a2 * t);
Since we are dividing by v(1), which is extremely small in your case:
>> v(:,1)
ans =
1.0e-03 *
0.005000000000000
0.004999999999140
0.004999999996559
0.004999999992258
0.004999999986237
0.004999999930331
0.004999999831432 ...
It might be treating it as a discontinuities in your ODE, and all solvers including ode45 and ode15 may be having problems.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by