0
$\begingroup$

I have problem with my code in Mathematica. I have introduced the set of coupled non-linear ODES. This is the resolution part:

(*Initial parameters*)
A = 0.5;
a = 0.9;
Ω = 0.24;

(*Initial conditions*)
υ0 = 0.22;
α0 = Pi;
ψ0 = Pi/2;
r0 = 20;
θ0 = Pi/8;
φ0 = 0;

Needs["DifferentialEquations`NDSolveProblems`"]; \
Needs["DifferentialEquations`NDSolveUtilities`"];
(*Systems to integrate*)

system = {x1'[t] == 
    Eq1[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
   x2'[t] == 
    Eq2[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
   x3'[t] == 
    Eq3[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
   x4'[t] == 
    Eq4[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
   x5'[t] == 
    Eq5[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]],
   x6'[t] == 
    Eq6[A, a, Ω, x1[t], x2[t], x3[t], x4[t], x5[t]], 
   x1[0] == υ0, x2[0] == α0, x3[0] == ψ0, 
   x4[0] == r0, x5[0] == θ0, x6[0] == φ0};

sol = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 14000}, 
   Method -> {"StiffnessSwitching", 
     Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 22,
    MaxSteps -> Infinity, PrecisionGoal -> 15, WorkingPrecision -> 22];

ParametricPlot3D[
 Evaluate[{x4[t]*Sin[x5[t]]*Cos[x6[t]], x4[t]*Sin[x5[t]]*Sin[x6[t]], 
    x4[t]*Cos[x5[t]]} /. sol], {t, 0, 14000}, PlotPoints -> 10000, 
 ColorFunction -> {Red}, ImageSize -> 500]

I receive the following error messages

NDSolve::precw: "The precision of the differential equation ({<<1>>}) is less than \ WorkingPrecision (22.`)"
NDSolve::ndsz: At t == 140.91450584595810589848638366914914657367`22., step size is \ effectively zero; singularity or stiff system suspected

Someone could suggest me how to improve my code? Thank you in advance.

$\endgroup$
13
  • 2
    $\begingroup$ Try to Rationailze your numbers, like A = Rationalize[0.5];. Does that hep? $\endgroup$
    – user21
    Sep 26 '19 at 10:44
  • $\begingroup$ Unfortunately not, I think it is more connected with the method to be used! :-( $\endgroup$
    – VDF
    Sep 26 '19 at 10:48
  • 3
    $\begingroup$ Please clarify. It you what @user21 suggests for all four parameters that are machine precision, the error (it's just a warning, actually) you cite goes away. Probably has to do with the missing code. Problems with code generally requires the code for the problem to be solved. You could try Rationalize[system] or Rationalize[system, 0]. $\endgroup$
    – Michael E2
    Sep 26 '19 at 10:57
  • $\begingroup$ It gives me still the same error, even if I add Rationalize $\endgroup$
    – VDF
    Sep 26 '19 at 12:22
  • $\begingroup$ I could not reproduce the test cases from the article arxiv.org/pdf/1901.03380.pdf . Perhaps there is an error in the equations. $\endgroup$ Sep 26 '19 at 22:15
2
$\begingroup$

I ran two tests to understand what freezes NDSolve. I used a simple code that implements Runge-Kutta 4 orders. As it turned out, NDSolve freezes at the moment when the particle descends from the regular trajectory. In the article on https://arxiv.org/abs/1901.03380v1 , the authors report "Therefore we adapted the highly-accurate core for the integration of photon trajectories used in LSDCode+ [45] to the case of massive particles. The code implements the Runge-Kutta method of the eighth order (the Dorman – Prince method) [46] with an adaptive step." I can advise you not to use NDSolve to solve the problem, but to develop your own code. Code for one test

(*Useful initial functions*)\[CapitalSigma][r_, a_, \[Theta]_] := 
  r^2 + (a*Cos[\[Theta]])^2;
\[CapitalDelta][r_, a_] := r^2 - 2 M*r + a^2;
\[Rho][r_, a_, \[Theta]_] := 
  r^2 + a^2 + 
   2 M*r*(a*Sin[\[Theta]])^2/\[CapitalSigma][r, a, \[Theta]];
\[Gamma][\[Upsilon]_] := 1/Sqrt[Abs[(1 - \[Upsilon]^2)]];

(*Useful further functions*)
N2[r_, a_, \[Theta]_] := \[CapitalDelta][r, a]/\[Rho][r, a, \[Theta]];
Nphi[r_, a_, \[Theta]_] := -2 M*a*
   r/(\[CapitalSigma][r, a, \[Theta]]*\[Rho][r, a, \[Theta]]);

(*Metric components*)
gtphi[r_, a_, \[Theta]_] := -4 M*a*r*
   Sin[\[Theta]]^2/\[CapitalSigma][r, a, \[Theta]];
gphiphi[r_, a_, \[Theta]_] := \[Rho][r, a, \[Theta]]*Sin[\[Theta]]^2;
grr[r_, a_, \[Theta]_] := \[CapitalSigma][r, 
    a, \[Theta]]/\[CapitalDelta][r, a];
gthth[r_, a_, \[Theta]_] := \[CapitalSigma][r, a, \[Theta]];
gtt[r_, a_, \[Theta]_] := -(1 - 2 M*r/\[CapitalSigma][r, a, \[Theta]]);

(*Kinematical quantities*)
AR[r_, a_, \[Theta]_] := (M/(\[Rho][r, a, \[Theta]]*
       Sqrt[\[CapitalSigma][r, a, \[Theta]]^5*\[CapitalDelta][r, 
          a]]))*(\[CapitalSigma][r, a, \[Theta]]^2*(r^2 - 
        a^2) + (a*
         Sin[\[Theta]])^2*(r^2*(3*r^2 - 4 M*r + 
           a^2) + (a*Cos[\[Theta]])^2*(r^2 - a^2)));
TR[r_, a_, \[Theta]_] := 
  a*M Sin[\[Theta]]*((r^2 + a^2)*(\[CapitalSigma][r, a, \[Theta]] - 
         2*r^2) - 
      2*r^2*\[CapitalSigma][r, a, \[Theta]])/(\[Rho][r, a, \[Theta]]*
      Sqrt[\[CapitalSigma][r, a, \[Theta]]^5]);
KR[r_, a_, \[Theta]_] := -Sqrt[\[CapitalDelta][r, 
       a]/\[CapitalSigma][r, 
        a, \[Theta]]^5]*(r*\[CapitalSigma][r, a, \[Theta]]^2 + 
      M (a*Sin[\[Theta]])^2*(\[CapitalSigma][r, a, \[Theta]] - 
         2*r^2))/(\[Rho][r, a, \[Theta]]);

AT[r_, a_, \[Theta]_] := -a^2*r M*
   Sin[2*\[Theta]]*(r^2 + a^2)/(\[Rho][r, a, \[Theta]]*
      Sqrt[\[CapitalSigma][r, a, \[Theta]]^5]);
TT[r_, a_, \[Theta]_] := 
  a^2*r M*Sin[2*\[Theta]]*Sin[\[Theta]]*
   Sqrt[\[CapitalDelta][r, a]]/(\[Rho][r, a, \[Theta]]*
      Sqrt[\[CapitalSigma][r, a, \[Theta]]^5]);
KT[r_, a_, \[Theta]_] := -Sin[
     2*\[Theta]]*((r^2 + 
         a^2)*(2*a^2*r M*
          Sin[\[Theta]]^2 + \[CapitalSigma][r, a, \[Theta]]^2) + 
      2*a^2*r M*\[CapitalSigma][r, a, \[Theta]]*
       Sin[\[Theta]]^2)/(2*\[Rho][r, a, \[Theta]]*
      Sqrt[\[CapitalSigma][r, a, \[Theta]]^5]*Sin[\[Theta]]^2);

(*Impact parameters and emission angles*)
RS = 5/2;(*radius of the emission source*)
b[a_, \[Theta]_, \[CapitalOmega]_] := -(gtphi[RS, a, \[Theta]] + 
     gphiphi[RS, a, \[Theta]]*\[CapitalOmega])/(gtt[RS, a, \[Theta]] +
     gtphi[RS, a, \[Theta]]*\[CapitalOmega]);
q[a_, \[Theta]_, \[CapitalOmega]_] := 
  If[b[a, \[Theta], \[CapitalOmega]] != 
    0, (b[a, \[Theta], \[CapitalOmega]]*Cot[\[Theta]])^2 - (a*
       Cos[\[Theta]])^2, -(a*Cos[\[Theta]])^2];
\[Beta][r_, a_, \[Theta]_, \[CapitalOmega]_] := 
  ArcCos[b[a, \[Theta], \[CapitalOmega]]*
    Sqrt[N2[r, 
       a, \[Theta]]]/(Sqrt[
        gphiphi[r, a, \[Theta]]]*(1 + 
         b[a, \[Theta], \[CapitalOmega]]*Nphi[r, a, \[Theta]]))];

(*Factor of the radiation field*)
Rrad[r_, a_, \[Theta]_, \[CapitalOmega]_] := (r^2 + a^2 - 
      a*b[a, \[Theta], \[CapitalOmega]])^2 - \[CapitalDelta][r, 
     a]*(q[a, \[Theta], \[CapitalOmega]] + (b[
          a, \[Theta], \[CapitalOmega]] - a)^2);
FACT[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := 
  A*(1 + b[a, \[Theta], \[CapitalOmega]]*Nphi[r, a, \[Theta]])^2/(N2[
       a, \[Theta], \[CapitalOmega]]*
      Sqrt[Rrad[r, a, \[Theta], \[CapitalOmega]]]);

(*Radiation field components*)
F1[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := 
  FACT[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi], 
    r, \[Theta]]*(1 - \[Upsilon]*Sin[\[Psi]]*
      Cos[\[Alpha] - \[Beta][r, 
         a, \[Theta], \[CapitalOmega]]])*(Sin[\[Psi]]*
      Cos[\[Alpha] - \[Beta][r, 
         a, \[Theta], \[CapitalOmega]]] - \[Upsilon]);

F2[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := 
  FACT[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi], 
    r, \[Theta]]*(1 - \[Upsilon]*Sin[\[Psi]]*
      Cos[\[Alpha] - \[Beta][r, a, \[Theta], \[CapitalOmega]]])*
   Cos[\[Psi]]*
   Cos[\[Alpha] - \[Beta][r, a, \[Theta], \[CapitalOmega]]]/\[Upsilon];

F3[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := 
  FACT[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi], 
    r, \[Theta]]*(1 - \[Upsilon]*Sin[\[Psi]]*
      Cos[\[Alpha] - \[Beta][r, a, \[Theta], \[CapitalOmega]]])*
   Sin[\[Alpha] - \[Beta][r, 
       a, \[Theta], \[CapitalOmega]]]/(\[Upsilon]*Sin[\[Psi]]);

(*EQUATIONS OF MOTION*)
Eq1[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := -1/\[Gamma][\[Upsilon]]*(Sin[\[Alpha]]*
       Sin[\[Psi]]*(AR[r, a, \[Theta]] + 
         2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]*TR[r, a, \[Theta]]) + 
      Cos[\[Psi]]*(AT[r, a, \[Theta]] + 
         2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]*TT[r, a, \[Theta]])) +
    F1[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi], 
    r, \[Theta]];

Eq2[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := \[Gamma][\[Upsilon]]/\[Upsilon]*(Sin[\[Psi]]*(AT[
          r, a, \[Theta]] + 
         2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]^2*TT[r, a, \[Theta]] +
          KT[r, a, \[Theta]]*\[Upsilon]^2*Cos[\[Alpha]]^2) - 
      Sin[\[Alpha]]*
       Cos[\[Psi]]*(AR[r, a, \[Theta]] + 
         2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]*TR[r, a, \[Theta]] + 
         KR[r, a, \[Theta]]*\[Upsilon]^2)) + 
   F2[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi], 
    r, \[Theta]];

Eq3[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := -\[Gamma][\[Upsilon]]*
    Cos[\[Alpha]]/(\[Upsilon]*Sin[\[Psi]])*(AR[r, a, \[Theta]] + 
      2*\[Upsilon]*Cos[\[Alpha]]*Sin[\[Psi]]*TR[r, a, \[Theta]] + 
      KR[r, a, \[Theta]]*\[Upsilon]^2 + 
      KT[r, a, \[Theta]]*\[Upsilon]^2*Cos[\[Psi]]^2*Sin[\[Alpha]]) + 
   F3[A, a, \[CapitalOmega], \[Upsilon], \[Alpha], \[Psi], 
    r, \[Theta]];

Eq4[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := \[Gamma][\[Upsilon]]*\[Upsilon]*Sin[\[Alpha]]*
   Sin[\[Psi]]/Sqrt[grr[r, a, \[Theta]]];

Eq5[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := \[Gamma][\[Upsilon]]*\[Upsilon]*
   Cos[\[Psi]]/Sqrt[gthth[r, a, \[Theta]]];

Eq6[A_, a_, \[CapitalOmega]_, \[Upsilon]_, \[Alpha]_, \[Psi]_, 
   r_, \[Theta]_] := \[Gamma][\[Upsilon]]*\[Upsilon]*Sin[\[Psi]]*
    Cos[\[Alpha]]/
     Sqrt[gphiphi[r, a, \[Theta]]] - \[Gamma][\[Upsilon]]*
    Nphi[r, a, \[Theta]]/Sqrt[N2[r, a, \[Theta]]];
(*Initial parameters*)
A = 85/100;(*These are a data of a simulation where A=0.85,a=0.07, \
\Omega=0.005,f=100. The initial conditions on the test particle \
are:r_0=15, \varphi_0=0, \theta_0=pi/4,\n\
u_0=0.13,\alpha_0=0,\psi_0=pi/2. I do not know how to copy my \
data.They are too many!*)
a = .5;
\[CapitalOmega] = .24;
M = 1;
(*Initial conditions*)
\[Upsilon]0 = .22;
\[Alpha]0 = Pi;
\[Psi]0 = Pi/2;
r0 = 20;
\[Theta]0 = Pi/8;
\[CurlyPhi]0 = 0;


system = {x1'[t] == 
    Eq1[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]], 
   x2'[t] == 
    Eq2[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]], 
   x3'[t] == 
    Eq3[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]], 
   x4'[t] == 
    Eq4[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]], 
   x5'[t] == 
    Eq5[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]], 
   x6'[t] == 
    Eq6[A, a, \[CapitalOmega], x1[t], x2[t], x3[t], x4[t], x5[t]], 
   x1[0] == \[Upsilon]0, x2[0] == \[Alpha]0, x3[0] == \[Psi]0, 
   x4[0] == r0, x5[0] == \[Theta]0, x6[0] == \[CurlyPhi]0};
sol = NDSolveValue[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 200}];

lst = Table[{x4[t]*Sin[x5[t]]*Cos[x6[t]], x4[t]*Sin[x5[t]]*Sin[x6[t]],
     x4[t]*Cos[x5[t]]}, {t, 0, 132, 1}];
(*RK4*) rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] := 
  Module[{table, ylist, step, k1, k2, k3, k4},
   step = N[(tfinal - tinit)/(nsteps)];
   ylist = valtinit;

   table = {ylist};
   Table[k1 = step*f /. MapThread[Rule, {variables, ylist}]; 
    k2 = step*f /. MapThread[Rule, {variables, k1/2 + ylist}];
    k3 = step*f /. MapThread[Rule, {variables, k2/2 + ylist}];
    k4 = step*f /. MapThread[Rule, {variables, k3 + ylist}];
    ylist += 1/6 (k1 + 2 (k2 + k3) + k4);
    AppendTo[table, ylist];
    ylist, nsteps];
   table];


funclist = {Eq1[A, a, \[CapitalOmega], x1, x2, x3, x4, x5], 
   Eq2[A, a, \[CapitalOmega], x1, x2, x3, x4, x5], 
   Eq3[A, a, \[CapitalOmega], x1, x2, x3, x4, x5], 
   Eq4[A, a, \[CapitalOmega], x1, x2, x3, x4, x5], 
   Eq5[A, a, \[CapitalOmega], x1, x2, x3, x4, x5], 
   Eq6[A, a, \[CapitalOmega], x1, x2, x3, x4, x5]};
initials = {\[Upsilon]0, \[Alpha]0, \[Psi]0, 
   r0, \[Theta]0, \[CurlyPhi]0};
variables = {x1, x2, x3, x4, x5, x6};
init = 0;


final = 2000; nstep = 20000; 
 sol4 = rk4[funclist, variables, initials, init, final, 
   nstep]; // AbsoluteTiming

st = N[(final - init)/(nstep)]; x1 = 
 Interpolation[Table[{i st, sol4[[i, 1]]}, {i, 1, nstep}]]; x2 = 
 Interpolation[Table[{i st, sol4[[i, 2]]}, {i, 1, nstep}]]; x3 = 
 Interpolation[Table[{i st, sol4[[i, 3]]}, {i, 1, nstep}]]; x4 = 
 Interpolation[Table[{i st, sol4[[i, 4]]}, {i, 1, nstep}]]; x5 = 
 Interpolation[Table[{i st, sol4[[i, 5]]}, {i, 1, nstep}]]; x6 = 
 Interpolation[Table[{i st, sol4[[i, 6]]}, {i, 1, nstep}]];

Show[ParametricPlot3D[
   Evaluate[{x4[t]*Sin[x5[t]]*Cos[x6[t]], x4[t]*Sin[x5[t]]*Sin[x6[t]],
      x4[t]*Cos[x5[t]]}], {t, 0, 2000}, PlotStyle -> Red, 
   PlotRange -> All], ListPointPlot3D[lst]] // Quiet

Figure 1 How to make NDSolve solve this problem? To do this, we define the classical Runge-Kutta method

ClassicalRungeKutta /: 
 NDSolve`InitializeMethod[ClassicalRungeKutta, __] := 
 ClassicalRungeKutta[]
ClassicalRungeKutta[___]["Step"[f_, t_, h_, y_, yp_]] := 
  Block[{deltay, k1, k2, k3, k4},
   k1 = yp;
   k2 = f[t + 1/2 h, y + 1/2 h k1];
   k3 = f[t + 1/2 h, y + 1/2 h k2];
   k4 = f[t + h, y + h k3];
   deltay = h (1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4);
   {h, deltay}
   ];

Then run three tests with a different combination of parameters

sol1 = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 2000}, 
    Method -> ClassicalRungeKutta, StartingStepSize -> 1/20, 
    MaxStepSize -> .05, MaxSteps -> 10^6]; // AbsoluteTiming
sol2 = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 2000}, 
    Method -> ClassicalRungeKutta, StartingStepSize -> 1/50, 
    MaxStepSize -> .02, MaxSteps -> 10^6]; // AbsoluteTiming
sol3 = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 2000}, 
    Method -> ClassicalRungeKutta, StartingStepSize -> 1/100, 
    MaxStepSize -> .01, MaxSteps -> 10^6, 
    WorkingPrecision -> 30]; // AbsoluteTiming

Comparing the three solutions, we see that the solution does not converge. This is the main reason NDSolve stops at t=132. using standard method. Figure 2

I found a combination of methods for rk8. Test case

A = 0.5;
a = 0.9;
\[CapitalOmega] = 0.24;

(*Initial conditions*)
\[Upsilon]0 = 0.22;
\[Alpha]0 = Pi;
\[Psi]0 = Pi/2;
r0 = 20;
\[Theta]0 = Pi/8;
\[CurlyPhi]0 = 0; tm = 14000;
sol1 = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, tm}, 
    Method -> {"FixedStep", "StepSize" -> .001, 
      Method -> {"ExplicitRungeKutta", 
        "DifferenceOrder" -> 8}}]; // AbsoluteTiming
(*Out[]= {380.243, Null}*)

ParametricPlot3D[
 Evaluate[{x4[t]*Sin[x5[t]]*Cos[x6[t]], x4[t]*Sin[x5[t]]*Sin[x6[t]], 
    x4[t]*Cos[x5[t]]} /. sol1], {t, 0, tm}, PlotStyle -> Red, 
 PlotRange -> All, AxesLabel -> {x, y, z}]

Figure 3

$\endgroup$
3
  • $\begingroup$ Thank you very much Alex for your help. So, you suggest me to develop my own code. I would have liked to do in Mathematica. Is it not possible to do at all? $\endgroup$
    – VDF
    Sep 27 '19 at 21:36
  • 1
    $\begingroup$ You saw in my code how I used rk4 to solve the problem. Similarly, you can build rk8. $\endgroup$ Sep 27 '19 at 23:53
  • $\begingroup$ The main problem of this code is Eq3. There is the first term $1 /\sin {\psi}$. When the system approaches $\psi = 0$ then NDSolve stops. I checked that rk8 does not go through $\psi = 0$. Interestingly, rk4 goes through $\psi= 0$. This is a classic! $\endgroup$ Sep 28 '19 at 13:00
0
$\begingroup$

I think it is the singularity of 1/x1[t] that stops NDSolve[], and that is an intrinsic difficulty.

In my experience, often "ImplicitRungeKutta" would be good enough for many stiff problems, as long as the ODE is solvable and analytical (mathematically infinite smooth).

sol = NDSolve[system, {x1, x2, x3, x4, x5, x6}, {t, 0, 200}, 
              Method -> "ImplicitRungeKutta"];

I copy-and-paste your code in both questions, then draw the curve of x1[t] with t until just before the equations blow up. It is hitting zero.

Plot[Evaluate[{x1[t]} /. sol], {t, 0, 140.9}]

enter image description here

Note that the Eq2 code reads

$$ \text{Eq2}[\text{A$\_$},\text{a$\_$},\Omega \_,\upsilon \_,\alpha \_,\psi \_,\text{r$\_$},\theta \_]\text{:=}\gamma [\upsilon ]/\upsilon *(\text{Sin}[\psi ]*(\text{AT}[r,a,\theta ]+2*\upsilon *\text{Cos}[\alpha ]*\text{Sin}[\psi ]{}^{\wedge}2*\text{TT}[r,a,\theta ]+\text{KT}[r,a,\theta ]*\upsilon {}^{\wedge}2*\text{Cos}[\alpha ]{}^{\wedge}2)-\text{Sin}[\alpha ]*\text{Cos}[\psi ]*(\text{AR}[r,a,\theta ]+2*\upsilon *\text{Cos}[\alpha ]*\text{Sin}[\psi ]*\text{TR}[r,a,\theta ]+\text{KR}[r,a,\theta ]*\upsilon {}^{\wedge}2))+\text{F2}[A,a,\Omega ,\upsilon ,\alpha ,\psi ,r,\theta ]; $$

Here is the problem: the term $\gamma [\upsilon ]/\upsilon$ ($\upsilon$ will be replaced by x1[t] above) hits the singularity point. That's why ODE solvers complain.

I don't have a good understanding of the ODE system here, but if the singularity is an illusion (say removable by reformulation or coordinate transformation), then there might be a chance to solve it.

$\endgroup$
5
  • 1
    $\begingroup$ The main problem of this code is Eq3. There is the first term 1/sinψ. When the system approaches ψ=0 then NDSolve stops. $\endgroup$ Sep 28 '19 at 13:07
  • $\begingroup$ @Alex Trounev Ah, I think you are right, it pass through zero several times already. $\endgroup$
    – Eddy Xiao
    Sep 28 '19 at 13:15
  • $\begingroup$ Thanks to all. The singularity you talk is removable, because it can be integrated. I should think how to recast it in a way that the apparent singularity is removed. Thanks for your suggestions. I will try to work on it $\endgroup$
    – VDF
    Sep 28 '19 at 17:35
  • $\begingroup$ I am still not able to solve this issue. The term 1/sin\psi is in all Eq. (3). Probably, I should put an "if", but how? $\endgroup$
    – VDF
    Sep 30 '19 at 9:15
  • $\begingroup$ I won't feel too bad if something like If[Abs[\Upsilon]<1/1000,...] is in the RHS of ODE, especially when you have to. I guess it will be helpful if you explain a bit about why the singularity is removable in the question. $\endgroup$
    – Eddy Xiao
    Oct 1 '19 at 10:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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