2
$\begingroup$

I want to solve this system of differential equations with Mathematica:

NDSolve[{y''[t] + 
3 y'[t] Sqrt[(y'[t])^2/6 + (1 - Exp[-y[t] Sqrt[2/3]])^2/4] + 
Sqrt[3/2] Exp[-y[t] Sqrt[2/3]] (1 - Exp[-y[t] Sqrt[2/3]]) == 0, 
q'[t]/q[t] == Sqrt[(y'[t])^2/6 + (1 - Exp[-y[t] Sqrt[2/3]])^2/4], 
m'[t] == 1/q[t], y'[0] == -0.008226306418212731, 
y[0] == 5.630991866033891, q[0] == 2.6791230019501455`*^-33, 
m[150] == 0}, {y, q, m}, {t, 0, 150}];

However, I am encountering with trouble as getting these messages:

Power::infy: Infinite expression 1/0. encountered.
Infinity::indet: Indeterminate expression (0. ComplexInfinity)/Sqrt[6] encountered.

Etc. And when trying to solve

sol = NDSolve[{y''[t] + 
3 y'[t] Sqrt[(y'[t])^2/6 + (1 - Exp[-y[t] Sqrt[2/3]])^2/4] + 
Sqrt[3/2] Exp[-y[t] Sqrt[2/3]] (1 - Exp[-y[t] Sqrt[2/3]]) == 0, 
q'[t]/q[t] == Sqrt[(y'[t])^2/6 + (1 - Exp[-y[t] Sqrt[2/3]])^2/4], 
y'[0] == -0.008226306418212731, 
y[0] == 5.630991866033891, q[0] == 2.6791230019501455`*^-33,}, {y, q}, {t, 0, 150}];
NDSolve[{m'[t] == 1/(q[t]/.sol[[1]]),m[150] == 0},m,{t,0,150}];

It says that the initial condition m[150]==1 is not compatible in dimensions with sol. So, how do I proceed?

$\endgroup$

2 Answers 2

7
$\begingroup$

Solve first for q[t],y[t] without using m[t] since m[t] is not coupled to the first two equations.

Then once you obtain q[t] solution, use that to solve for m[t] using the third ode on its own.

ode1 = y''[t] + 
   3  y'[t]  Sqrt[(y'[t])^2/6 + (1 - Exp[-y[t]  Sqrt[2/3]])^2/4] + 
   Sqrt[3/2]  Exp[-y[t]  Sqrt[2/3]]  (1 - Exp[-y[t]  Sqrt[2/3]]) == 0

ode2 = q'[t]/q[t] == 
  Sqrt[(y'[t])^2/6 + (1 - Exp[-y[t]  Sqrt[2/3]])^2/4]

ic = {y'[0] == -0.008226306418212731, y[0] == 5.630991866033891, 
  q[0] == 2.6791230019501455`*^-33}

sol = NDSolve[{ode1, ode2, ic}, {y, q}, {t, 0, 150}]

Now solve for m[t]

ode3 = m'[t] == 1/First@Evaluate[q[t] /. sol]
sol = NDSolve[{ode3, m[150] == 0}, m, {t, 0, 150}]

ps. stackexchange is broke now, can't paste images and SE uploader is also not working. So can't show output as images or show plots any more. Pasting images gives error that image is too small.

$\endgroup$
4
$\begingroup$

We can solve equations in two steps as Nasser mentioned, using NItegrate as follows

 sol = NDSolve[{y''[t] + 
      3 y'[t] Sqrt[(y'[t])^2/6 + (1 - Exp[-y[t] Sqrt[2/3]])^2/4] + 
      Sqrt[3/2] Exp[-y[t] Sqrt[2/3]] (1 - Exp[-y[t] Sqrt[2/3]]) == 0, 
    q'[t] == 
     q[t] (Sqrt[(y'[t])^2/6 + (1 - Exp[-y[t] Sqrt[2/3]])^2/4]), 
    y'[0] == -0.008226306418212731, y[0] == 5.630991866033891, 
    q[0] == 2.6791230019501455`*^-33}, {y, q}, {t, 0, 150}];

Visualization

{Plot[Evaluate[y[t] /. sol[[1]]], {t, 0, 150}, 
  AxesLabel -> {"t", "y"}], 
 Plot[Evaluate[q[t] /. sol[[1]]], {t, 0, 150}, PlotRange -> All, 
  AxesLabel -> {"t", "q"}]}

Figure 1 Then we define

m0 = NIntegrate[1/q[t] /. sol[[1]], {t, 0, 150}, 
  Method -> "LocalAdaptive"]; 
m[T_?NumericQ] := 
 NIntegrate[1/q[t] /. sol[[1]], {t, 0, T}, 
   Method -> "LocalAdaptive"] - m0;

Visualization

Plot[Evaluate[m[t]], {t, 0, 150}, PlotRange -> All]

Figure 2

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.