numerical solution of ode singularity

469 Views Asked by At

I have to solve a nonlinear ode problem: \begin{align} \dot\gamma&=-2(D+0.5A+0.5A\cos\gamma)\sin k \\ -\dot k&=(E-B-B\cos\gamma)+(2D+A+A\cos\gamma)\frac{\cos\gamma}{\sin\gamma}\cos k \end{align} with the initial condition $\gamma=\pi$ and $k=-\frac{\pi}{2}$.
My problem is that right hand side of the second equation is infinity for the initial condition, so I cannot use an Euler scheme, for example.

Is there a way to solve it? Can anyone recommend a good textbook on this kind of problem? Thanks a lot!

1

There are 1 best solutions below

0
On

Playing around with the equations a bit, you can observe that, introducing \begin{equation} X = \cos k\,\sin \gamma,\quad Y = \sin k\,\sin \gamma,\quad Z = \cos \gamma, \tag{1} \end{equation} we obtain the following system for $(X,Y,Z)$: \begin{align} \dot{X} &= \phantom{-}\left(E - B(1+Z)\right)\,Y,\\ \dot{Y} &= -\left(E - B(1+Z)\right)\,X - \left(2 D+A(1+Z)\right)\,Z,\tag{2}\\ \dot{Z} &= \phantom{-}\left(2 D+A(1+Z)\right)\,Y. \end{align} You can easily deduce that $\frac{\text{d}}{\text{d} t}\left(X^2 + Y^2 + Z^2 \right) = 0$, which makes perfect sense, since from $(1)$ follows that \begin{equation} X^2 + Y^2 + Z^2 = 1. \end{equation} In other words, system $(2)$ is naturally confined to the sphere in $(X,Y,Z)$-space.

Admittedly, it looks a bit strange to increase the number of model variables, but this serves a purpose. The initial condition $k(0) = -\frac{\pi}{2}, \gamma(0) = \pi$ now translates to \begin{equation} (X(0),Y(0),Z(0)) = (0,0,-1). \end{equation} In other words, the initial condition is at the 'south pole' of the sphere. This immediately explains why this initial condition gives infinities when used in the original system: the south pole is a singular point in the spherical coordinates $(1)$, see also MathWorld's lemma on spherical coordinates.

This suggests a solution to your problem. Note that the system $(2)$ is perfectly well-behaved at $(X,Y,Z) = (0,0,-1)$; the right hand side at this point evaluates to $(0,2 D,0)$. As a matter of fact, system $(2)$ behaves nicely everywhere. So, as it turns out, we just made an unfortunate choice of spherical coordinates to study solutions to your original system with this particular initial condition. Let's change our spherical coordinates then! A simple solution is to rotate the singular point, the south pole, away over a $90^\circ$ angle. We can do this by introducing \begin{equation} X = \cos \alpha\, \sin \beta,\quad Y = \cos \beta,\quad Z = \sin \alpha \sin \beta, \end{equation} which shifts the 'poles' of the coordinate system to $(0,\pm 1,0)$. The original system transforms to \begin{align} \dot{\alpha} &= \left(A \cos\alpha + B \sin \alpha\right) \sin\alpha\,\cos\beta + \left((A+2 D)\cos \alpha + (B-E)\sin \alpha\right) \frac{\cos \beta}{\sin\beta},\\ \dot{\beta} &= \left(2 D + A(1+\sin\alpha\,\sin\beta)\right)\sin\alpha + \left(E - B(1+\sin\alpha\sin\beta)\right)\cos\alpha. \end{align} The initial condition $(X,Y,Z) = (0,0,-1)$ is now uniquely transformed in the initial condition $\alpha = -\frac{\pi}{2}, \beta = \frac{\pi}{2}$. You should have no problems treating this system numerically.