12
$\begingroup$

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!

$\endgroup$
  • 2
    $\begingroup$ 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... $\endgroup$ – J. M.'s discontentment Aug 15 '15 at 13:01
18
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

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).

enter image description here
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}
 ]
| improve this answer | |
$\endgroup$
  • $\begingroup$ Very instructive and elegant. $\endgroup$ – chris Aug 16 '15 at 11:40
  • $\begingroup$ 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? $\endgroup$ – user1337 Nov 2 '17 at 10:13
8
$\begingroup$

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

enter image description here

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

enter image description here

| improve this answer | |
$\endgroup$

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.