TransWikia.com

Avoiding/skipping singularity error in ParametricNDSolve

Mathematica Asked on June 27, 2021

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.

One Answer

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.

Answered by Ulrich Neumann on June 27, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP