1
$\begingroup$

I've looked up a few examples of similar errors for ODE's with singular points, but have yet to figure out a way to solve this problem for my case. A simple patch up where the singular point is 'skipped' would be sufficient for me, but I would really appreciate some help.

Edit: my point being that ParametricNDSolve is giving the ParametricNDSolve::ndsz: At t == 0.13801826916665527`, step size is effectively zero; singularity or stiff system suspected. error

So I have a velocity equation V(x,t) for a particular as a function of (x,t) and some other parameters:

Edit: in my original unedited post, I think the code was copying across incorrectly. The following should show the issues with the ParametricNDSolve hitting the singularities.

ClearAll["Global`*"]

v = 0.3;
k0r = k0;
k0l = k0;
\[Sigma]r = \[Sigma];
\[Sigma]l = \[Sigma];
A[x_, t_, k0_, \[Sigma]_, \[Alpha]_] = 
  Sqrt[2/\[Pi]] \[Sigma]r Sqrt[(1 - v)/(
    1 + v)] \[Alpha] E^(-2 ((1 - v)/(1 + v)) (t - x)^2 \[Sigma]r^2) - 
   Sqrt[2/\[Pi]] \[Sigma]l Sqrt[(1 + v)/(
    1 - v)] (1 - \[Alpha]) E^(-2 ((1 + v)/(
      1 - v)) (t + x)^2 \[Sigma]l^2) + 
   Sqrt[(1 + v)/(
    1 - v)] (-Sqrt[(2/\[Pi])] Sqrt[\[Sigma]r \[Sigma]l] Sqrt[k0l/k0r]
       Sqrt[\[Alpha] (1 - \[Alpha])]
       E^(-((1 + v)/(1 - v)) (t + x)^2 \[Sigma]l^2 - ((1 - v)/(
         1 + v)) (t - x)^2 \[Sigma]r^2) (Cos[
         k0r (-Sqrt[((1 - v)/(1 + v))] (t - x)) + 
          k0l Sqrt[(1 + v)/(1 - v)] (t + x)] - 
        2*\[Sigma]l^2/k0l Sqrt[(1 + v)/(
         1 - v)] (t + x) Sin[
          k0r (-Sqrt[((1 - v)/(1 + v))] (t - x)) + 
           k0l Sqrt[(1 + v)/(1 - v)] (t + x)])) + 
   Sqrt[(1 - v)/(
    1 + v)] (Sqrt[2/\[Pi]] Sqrt[\[Sigma]r \[Sigma]l] Sqrt[k0r/k0l]
       Sqrt[\[Alpha] (1 - \[Alpha])]
       E^(-((1 + v)/(1 - v)) (t + x)^2 \[Sigma]l^2 - ((1 - v)/(
         1 + v)) (t - x)^2 \[Sigma]r^2) (Cos[
         k0r (-Sqrt[((1 - v)/(1 + v))] (t - x)) + 
          k0l (Sqrt[(1 + v)/(1 - v)] (t + x))] + 
        2*\[Sigma]r^2/k0r Sqrt[(1 - v)/(
         1 + v)] (t - x) Sin[
          k0r (-Sqrt[((1 - v)/(1 + v))] (t - x)) + 
           k0l Sqrt[(1 + v)/(1 - v)] (t + x)]));

B[x_, t_, k0_, \[Sigma]_, \[Alpha]_] = 
  Sqrt[2/\[Pi]] \[Sigma]r Sqrt[(1 - v)/(
    1 + v)] \[Alpha] E^(-2 ((1 - v)/(1 + v)) (t - x)^2 \[Sigma]r^2) + 
   Sqrt[2/\[Pi]] \[Sigma]l Sqrt[(1 + v)/(
    1 - v)] (1 - \[Alpha]) E^(-2 ((1 + v)/(
      1 - v)) (t + x)^2 \[Sigma]l^2) + 
   Sqrt[(1 + v)/(
    1 - v)] (Sqrt[2/\[Pi]] Sqrt[\[Sigma]r \[Sigma]l] Sqrt[k0l/k0r]
       Sqrt[\[Alpha] (1 - \[Alpha])]
       E^(-((1 + v)/(1 - v)) (t + x)^2 \[Sigma]l^2 - ((1 - v)/(
         1 + v)) (t - x)^2 \[Sigma]r^2) (Cos[
         k0r (-Sqrt[((1 - v)/(1 + v))] (t - x)) + 
          k0l Sqrt[(1 + v)/(1 - v)] (t + x)] - 
        2*\[Sigma]l^2/k0l Sqrt[(1 + v)/(
         1 - v)] (t + x) Sin[
          k0r (-Sqrt[((1 - v)/(1 + v))] (t - x)) + 
           k0l Sqrt[(1 + v)/(1 - v)] (t + x)])) + 
   Sqrt[(1 - v)/(
    1 + v)] (Sqrt[2/\[Pi]] Sqrt[\[Sigma]r \[Sigma]l] Sqrt[k0r/k0l]
       Sqrt[\[Alpha] (1 - \[Alpha])]
       E^(-((1 + v)/(1 - v)) (t + x)^2 \[Sigma]l^2 - ((1 - v)/(
         1 + v)) (t - x)^2 \[Sigma]r^2) (Cos[
         k0r (-Sqrt[((1 - v)/(1 + v))] (t - x)) + 
          k0l (Sqrt[(1 + v)/(1 - v)] (t + x))] + 
        2*\[Sigma]r^2/k0r Sqrt[(1 - v)/(
         1 + v)] (t - x) Sin[
          k0r (-Sqrt[((1 - v)/(1 + v))] (t - x)) + 
           k0l Sqrt[(1 + v)/(1 - v)] (t + x)]));

which I am using ParametricNDSolve to solve:

params = {x, t, k0, \[Sigma] , \[Alpha]};

V[x_, t_, 
   k0_, \[Sigma]_, \[Alpha]_] = (A[x, t, k0, \[Sigma], \[Alpha]])/(B[
     x, t, k0, \[Sigma], \[Alpha]]);

Startpoints[k0i_, \[Sigma]i_, vi_, \[Alpha]i_, fr_, to_, st_, 
  nu_] := (samps = 
   Table[B @@ (params /. {k0 -> k0i, \[Sigma] -> \[Sigma]i, 
         v -> vi, \[Alpha] -> \[Alpha]i} /. t -> -2), {x, fr, to, st}];
  tab = Table[Plus @@ #[[1 ;; n]]*st, {n, 1, Length[#]}] &[samps];
  Flatten[
     Table[FirstPosition[tab, x_ /; x > n, {Length[tab]}] - 1, {n, 0, 
       1, 1/(nu - 1)}]]*st + fr)


subs = {k0 -> 5, \[Sigma] -> 1, v -> 0, \[Alpha] -> Sqrt[1]/2};
de = {xtraj'[t] == V @@ (params /. subs /. {x -> xtraj[t]}), 
   xtraj[-2] == x0};
parasols = ParametricNDSolve[de, xtraj[t], {t, -2, 2}, x0];
startpoints = 
  Startpoints @@ ({k0, \[Sigma], v, \[Alpha], -3, 3, 0.01, 50} /. 
     subs);
trajectories = 
 Table[{xtraj[t][x0], t} /. parasols, {x0, 
   startpoints}]; ParametricPlot[trajectories, {t, -2, 2}, 
 AspectRatio -> 1/2, PlotStyle -> Red]

There are singular points in the V(x,t) equation, most obviously when the denominator equals zero.

Is there a way for me to tell Mathematica to 'skip' the singular points? At the singularity itself I expect the trajectory to go horizontal (the velocity being infinite for that point) before going backwards for a while, turning around and then hitting another singular point and continuing on. If one runs the code above, you can see that it Mathematica doesn't want to plot the initial conditions for which there are singular points.

$\endgroup$
1
$\begingroup$

StreamPlot shows the trajectories without solving the ode:

StreamPlot[{1, V @@ (params /. subs)}, {t, -2, 2}, {x , -3, 3},FrameLabel -> { t, x}]

enter image description here

No singular points to be recognized.

$\endgroup$
7
  • $\begingroup$ Hi Ulrich, doesn't StreamPlot just plot the velocity field rather than the trajectory x(t) parametrically? Also, it's a bit hard to tell whether the field shown above has singularities, but I have included a graph showing that the denominator of V(x,t) is definitely zero at points. Ideally I'd like to solve the ODE while dealing with the singular points. $\endgroup$
    – j.foobles
    May 3 at 11:41
  • $\begingroup$ @j.foobles You might plot your "trajectories" together with the stream plot, fits quite well! $\endgroup$ May 3 at 12:33
  • $\begingroup$ @j.foobles If I plot B[0,t,5,1,1/2] in the range -2<t<2 I see a function without zeros! $\endgroup$ May 3 at 12:39
  • $\begingroup$ Hi @Ulrich Neumann, I am a bit confused - perhaps some of the code copied across incorrectly from Mathematica. Nevertheless I am convinced that the denominator function has some zeroes, so I've attached a .nb file here: 1drv.ms/u/s!Ak_7kWMZ0x15hJpuG74r3VEUFQIkTQ?e=dbay6L $\endgroup$
    – j.foobles
    May 3 at 12:53
  • $\begingroup$ @j.foobles You evaluate the same streamplot as in my answer? $\endgroup$ May 3 at 12:57

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.