Numerical integration algorithm for this set of ODEs

Question

I have to solve numerically the following initial value problem, which physical origin is discussed in this article by G.I. Taylor (1941). Here $\eta$ is the independent variable and ranges beetween $0$ and $1$, $\gamma$ is fixed and $=1.4$ in my simulation:

$$f'=\dfrac{f[-3\eta + (3+\gamma/2)\phi - 2\gamma \dfrac{\phi ^2}{\eta}]}{(\eta-\phi)^2-f/\psi},$$

$$\phi'=\dfrac{\frac{1}{\gamma} f'}{\psi(\eta-\phi)}-\dfrac{3}{2}\dfrac{\phi}{\eta-\phi},$$

$$\psi ' = \dfrac{\phi ' + 2 \frac{\phi}{\eta}}{\eta-\phi}\psi.$$

Note that the equations are in normal form. The initial conditions are: $$\psi (1) = \dfrac{\gamma + 1}{\gamma -1},\quad f(1)=\dfrac{2\gamma}{\gamma + 1},\quad \phi (1)=\dfrac{2}{\gamma + 1}.$$

The result I'm expecting is more or less this, where density=$\psi$, velocity=$\phi$, pressure=$f$:

As you can see, $\phi$ becomes soon linear, with a slope that can be seen to be $\frac{1}{\gamma}$, while $\psi$ and $f$ respectively approach the costant values $0$ and (according to Taylor) $0.436$.

The integration scheme I'm using right now is the midpoint method: $$y_{n+1}=y_n + hF(x_n+\frac{h}{2},y_n+\frac{h}{2}F(x_n,y_n)).$$

Everything seems to work quite well until $\eta \approx 0.5$, $\phi'$ seems to approach $0.714=1/\gamma$ (however $f$ at that point is about $0.427$, which is a little lower than its limit and this should not occur). About at that point, $\phi'$ starts again to rise and for $\eta \to 0$ $\phi',\phi,\psi'$ and $\psi$ reach astronomical values. This is a plot of $\phi$ vs. $\eta$; here the ihtegration was done with a step $h=0.0001$:

I think that the main problem lies in the approximation errors due to the little numbers that enter the ratios I have to calculate in order to compute $f',\phi ',\psi '$ (that are computed as in the equations above in my program). If it can be of any relevance, I'm writing the program in C.

Could you suggest any improvement in the algorithm (or maybe a more robust one) and/or the possible problems which may arise in this computation?

---

I have tried the substitution $\xi = \phi - \eta/\gamma$ as suggested by Kirill and effectively the system seems to be more stable. However I'm still getting an explosion near $0$:

Answer 1

It is entirely possible, although I cannot tell by looking at it, that your system becomes stiff near $\eta=0$. I say this because Mathematica's NDSolve integrates this up to $\eta=3.4\times 10^{-9}$, then complains about a singularity or stiffness.

The midpoint method is not A-stable, so that may be the source of inaccuracy: it is unsuitable for stiff problems. The ODE has some kind of singularity near $0$, so maybe the solution blows up because it hits the singularity due to inaccuracy before the method's instability becomes apparent. What do you get if you use a larger time step, like $0.01$ or $0.001$?

Try an A-stable method, like the trapezoid rule: $$ y_{n+1} = y_n + \frac{h}{2}(F(x_n,y_n) + F(x_{n+1}, y_{n+1})). $$ You'd need to solve the implicit system of nonlinear equations here.

Another approach, since you know the asymptotic behaviour of the ODE solution, is to write the unknown functions as perturbations around the known asymptotic behaviour (e.g., $\frac1\gamma\eta+\eta^2\hat\phi(\eta))$, etc., as in a Taylor series), which might make the resulting ODE in $\hat\phi$ easier to solve by (maybe) avoiding the numerical errors in the denominator.

Answer 2

Instead of trying to shoehorn this into an ODE solver you write yourself (that probably uses a simplistic time step length selection method, low order, etc), use one of the existing, advanced and very accurate solvers for ODEs that you can find in any software package -- e.g., dsolve in Maple, NDSolve in Mathematica, or ode45 in Matlab. If necessary, play with the parameters these functions have and that you can glean from their respective help pages.

The point being that these are fantastic methods that give you high accuracy you could never get from a method you implement yourself in a finite amount of time. There's no point wasting your time reinventing something others have already done.

Answer 3

My suggestion is to analytically detect the singularities and then to switch form explicit solver to implicit one around the singularity. From this equation $ \psi ' = \dfrac{\phi ' + 2 \frac{\phi}{\eta}}{\eta-\phi}\psi. $ it follows that the singularity is around $\eta-\phi$. Therefore, you should require that $ \phi ' = - 2 \frac{\phi}{\eta} $ with a solution $\phi=\frac{c}{{\eta}^{2}}$. Therefore $ \eta = \frac{c}{{\eta}^{2}} $ and the critical values is $\eta^*= ^3\sqrt{c}$.

The second equation $ \phi'=\dfrac{\frac{1}{\gamma} f'}{\psi(\eta-\phi)}-\dfrac{3}{2}\dfrac{\phi}{\eta-\phi} $ will be then $ f' = \dfrac{3 \gamma \psi}{2} \phi $

The first equation $ f'=\dfrac{f[-3\eta + (3+\gamma/2)\phi - 2\gamma \dfrac{\phi ^2}{\eta}]}{(\eta-\phi)^2-f/\psi} $ will become $ f'=\psi(-3\eta + (3+\gamma/2)\phi - 2\gamma \dfrac{\phi ^2}{\eta}) $