# Extending NDSolve beyond a singularity

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},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!

• Interesting question! A very different approach will be needed, since InterpolatingFunction[] assumes that whatever it's approximating is continuous (thus, pole-free) throughout its domain... – J. M.'s ennui Aug 15 '15 at 13:01

## 2 Answers

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}
]

• Very instructive and elegant. – chris Aug 16 '15 at 11:40
• Thanks for the lovely method Michael. Although it has been years since the question was asked, I've been recently trying to extend this method to ODE systems and/or higher order ODEs, such as the Weierstrass equation $y''=6y^2-1/2g_2$. However, I couldn't do it by embedding $(y,y')$ in either $\mathbb{RP}^1 \times \mathbb{RP}^1$ nor $\mathbb{RP}^2$. Any tips on how to extend your method to such systems? – user1337 Nov 2 '17 at 10:13

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,
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 == 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 == 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}}] 