Extending NDSolve beyond a singularity

Question

The $\tan$ function satisfies the following IVP:

$$y'=1+y^2 ,\quad y(0)=0 $$

and has simple poles at the points $x=\pi/2+ \pi n$ for integer $n$.

When trying to get $\tan$ via numerical integration, the command

NDSolve[{y'[x]==y[x]^2+1,y[0]==0},y[x],{x,-10,10}]

gives a solution which is defined only for $x \in(- \pi/2,\pi/2)$. Is there a way to extend the solution beyond the poles $x= \pm \pi/2$? What about singularities in the general case?

Thank you!

Answer 1

We can treat the variable $y$ as an element $[y_1 \colon y_2]$ of the projective line. In code, this means replacing y[x] by y1[x]/y2[x]. For an IVP $y' = f(x, y), \ y(x_0) = y_0$, we translate the initial condition as $y_1(x_0) = y_0, \ y_2(x_0) = 1$. Since the substitution yields an equation in two variables $y_1$, $y_2$, $$y_1'y_2-y_1y_2'=y_2^2\;f(x,\,y_1/y_2)\,,$$ we need another equation. So to get a unique solution, we impose the condition in code as

y1[x]^2 + y2[x]^2 == y0^2 + 1

Since this condition is satisfied initially, NDSolve will use it in conjunction with the ODE to determine y1[x] and y2[x] at each step. We can use this condition as it is and solve the system as a differential-algebraic equation (DAE); or we can differentiate it and solve the system as an ODE. The important difference is that the methods and precision available for DAEs are limited.

eqn = y'[x] == 1 + y[x]^2;
blowup = {y -> (y1[#]/y2[#] &)};

newfn = eqn /. blowup /. Equal -> Subtract // Together // Numerator;

newDAE = {newfn == 0, y1[x]^2 + y2[x]^2 == y0^2 + 1};
newODE = {newfn == 0, D[y1[x]^2 + y2[x]^2 == y0^2 + 1, x]};

Block[{x0 = 0, y0 = 0},
  sol = NDSolve[{newDAE, y1[x0] == y0, y2[x0] == 1}, {y1, y2}, {x, 0, 10}]
  ];

Block[{x0 = 0, y0 = 0},
  sol = NDSolve[{newODE, y1[x0] == y0, y2[x0] == 1}, {y1, y2}, {x, 0, 10}]
  ];

Both solutions yield the same plots:

Plot[y[x] /. blowup /. First@sol // Evaluate, {x, 0, 10}]

Compare by overlaying the graph of tangent, the exact solution in this example:

Plot[{y[x] /. blowup /. First@sol, Tan[x]} // Evaluate, {x, 0, 10}]

Update: Another view of what is happening.

A standard model of the projective line $[y_1\colon y_2]$ is a unit-diameter circle tangent to an axis. The corresponding affine line $y$ is given by $y_2 = 1$. Here we project the solution in terms of {y1[x], y2[x]} onto the desired solution y[x] (for x running from 0 to 10).

The projection from the circle model of the projective line onto the affine line y2 == 1. (The cylinder is the product of the interval 0 <= x <= 10 and the projective line or circle.)

cplot2 = ContourPlot3D[y1^2 + y2^2 == y2,
   {x, 0, 10}, {y1, -1.05, 1.05}, {y2, -0.05, 1.05},
   ContourStyle -> Opacity[0.3], Mesh -> None];
base = Show[
   ParametricPlot3D[
    Evaluate[{x, ( y1[x] y2[x])/(y1[x]^2 + y2[x]^2), y2[x]^2/(
       y1[x]^2 + y2[x]^2)} /. First@sol], {x, 0, 10}],
   ParametricPlot3D[
    Evaluate[{x, y[x] /. blowup, 1} /. First@sol], {x, 0, 10}, 
    PlotStyle -> ColorData[97, 3], Exclusions -> Cos[x] == 0],
   (*cplot1,*)cplot2,
   PlotRange -> {{0, 10}, {-2, 2}, {-0.1, 3.05}}, 
   AxesLabel -> {x, y1, y2}];
(* * * * *)
Manipulate[
 Show[
  base,
  Graphics3D[{
    Gray,
    Table[InfiniteLine[{{0, y, 1}, {10, y, 1}}], {y, -2, 2}],
    Table[
     InfiniteLine[{{x0, -1, 1}, {x0, 1, 1}}], {x0, 0, 10, Pi/2}],
    Red, Thickness[Medium],
    Line[{{0, 0, 0}, {10, 0, 0}}],
    InfiniteLine[{{x, 0, 0}, {x, y[x] /. blowup /. First@sol, 1}}],
    PointSize[Large], 
    Point[{{x, 0, 0}, {x, y[x], 1}, {x, ( y1[x] y2[x])/(
         y1[x]^2 + y2[x]^2), y2[x]^2/(y1[x]^2 + y2[x]^2)}} /. 
       blowup /. First@sol]
    }]
  ],
 {x, 0, 10}
 ]

Answer 2

You can use WhenEvent[ ] function, maybe it's not the perfect (not very accurate) solution but it works.

sol = With[{eps = 10^5}, 
First@NDSolve[{y'[x] == y[x]^2 + 1, y[0] == 0, 
WhenEvent[{y[x] == eps, y[x] == -eps}, y[x] -> -y[x]]}, 
y, {x, -10, 10}]];

Plot[{y[x] /. sol, Tan[x]}, {x, -6, 6}, 
PlotLegends -> {"NDSolve", "Tan[x]"}, 
PlotStyle -> {Thin, {Thick, Dashed}}]

Another example:

 sol2 = With[{eps = 10^5}, 
 First@NDSolve[{y'[x] == y[x]^2 + x, y[0] == 1, 
 WhenEvent[{y[x] == eps, y[x] == -eps}, y[x] -> -Abs[y[x]]]}, 
 y, {x, -10, 10}]];
 sol3 = First@DSolve[{y'[x] == y[x]^2 + x, y[0] == 1}, y[x], x];

 Plot[{y[x] /. sol2, y[x] /. sol3}, {x, 0, 10}, 
 PlotLegends -> {"NDSolve", "DSolve"}, 
 PlotStyle -> {Thin, {Thick, Dashed}}, 
 PlotRange -> {{0, 10}, {-10, 10}}]