How to solve a differential equation (with an indeterminate form at a point) using Runge-Kutta method?

Question

I am trying to plot the integral curves corresponding to the following differential equation:

where 'u' and 'a' are respectively the velocity of fluid and sound speed respectively.

The sonic point (r_s) of the flow is defined as the point where the flow velocity equals the sound speed obtained by setting both the numerator and denominator to zero. In the plot below, the sonic point is the point where the two transonic curves intersect.

Now I need to plot the Mach number (M=u/a) as a function of the radial coordinate to obtain the following plot:

I tried the following:

I tried to solve the differential equation using Runge-Kutta method. Since du/dr is 0/0 form at the sonic point (r_s), we need to use the L'Hospital's rule. But I could not understand how to deal with this.

Can anyone please suggest me how to proceed?

Answer 1

Your differential equation has the implicit solutions

$$ \frac{1}{r}+2\,{a}^{2}\ln \left( r \right) - \frac{ u ^{2}}{2}+{a}^{2}\ln \left( u \right) = c $$

Presumably we want $a, r, u > 0$. The minimum value of $1/r + 2 a^2 \ln(r)$ is $2 a^2 (1 - \ln(2 a^2))$, occurring at $r = 1/(2 a^2)$. The maximum value of $-u^2/2 + a^2 \ln(u)$ is $-a^2/2 + a^2 \ln(a)$, occurring at $u=a$. The critical value of the constant $c$ is $$c_0 = \frac{3 a^2}{2} - a^2 \ln(4 a^3)$$ If $c > c_0$, there is no solution for $u$ at $r=1/(2a^2)$.

I think there may be great difficulty trying to use Runge-Kutta for the solutions at $c = c_0$, because

1. At $r=1/(2a^2), u=a$, not only is the equation singular, but uniqueness breaks down as two solutions cross.
2. Arbitrarily close to this point, you get solutions (for $c > c_0$) that cease to exist. Any numerical method is likely to stray into this region.

Here is a plot of solutions for $c=c_0$ and some on each side, obtained using Maple in the case $a=1$.

What you might do is use the theory to approximate the solutions in the neighbourhood of the crossing point $(r=1/(2a^2), u=a)$, and only use Runge-Kutta (forward or backward) when safely away from this point.