ODE45: doubts about the result. Correct or not?

1 visualización (últimos 30 días)
Marco Sammito
Marco Sammito el 28 de Nov. de 2016
Comentada: bio lim el 28 de Nov. de 2016
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 comentario
bio lim
bio lim el 28 de Nov. de 2016
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.

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by