# Numerical integration algorithm for this set of ODEs

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$:

• Why is the first plot over $x\in[0,2]$, not $\eta\in[0,1]$? Can you post a plot that directly shows the solution blow up? Feb 8 '15 at 12:31
• @Kirill , thank you for your help, I will update soon following your suggestions. The plot is over $x\in (0,2)$ probably to evidence the discontinuity through the wave front - outside $(0,1)$ the solution is constant as shown in those plots, however I'm solving the system in $(0,1)$. Moreover, I'm not 100% sure about the definitions of the variables in the plot (I suppose that "pressure" is "differential pressure"), however they give a qualitative idea of what the solution should be. Feb 8 '15 at 15:38
• Also, in my answer I assumed you specifically wanted to implement your own method. If you are okay with using an existing implementation, Wolfgang Bangerth's advice is better. Feb 9 '15 at 0:13

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.

• Hi @Kirill, I tried to work with $\phi = \eta/\gamma + \xi$, and the system for $f,\xi,\psi$ seems to be more stable, however the explosion still occurs (I've posted a screenshot). About the trapezoid rule, how could I implement it in my code? (I mean, how could I solve the implicit system?) Feb 9 '15 at 14:25
• @pppqqq I specifically meant $\phi=\eta/\gamma+\eta^2\xi$ (maybe $\phi=\eta\xi$ would also work), so that small absolute errors in $\xi$ don't translate into large relative errors in $\phi$, unlike in your substitution. You can try solving it with a plain iterative method, or Newton's method. If you can't calculate the Jacobian explicitly, using finite differences should also work. Feb 9 '15 at 22:21
• I'm working right now at the code for $\phi = \eta /\gamma +\alpha * \eta ^ n *{\xi}$, with $n$ variable. Now I have the practical problem of writing the equations for $f'$ etc. If I have significant updates I'll post. By the way, your answer has been very helpful in understanding what's going on, thank you. Feb 10 '15 at 10:03
• @pppqqq I recommend using a CAS. Feb 10 '15 at 10:11

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.

• Hi @Wolfgang, totally aware of that, but my task is to write my own program [I have to do this for an exam (and it is fun!)]. Feb 9 '15 at 9:16

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})$

• It would be an improvement if you sketched how "to analytically detect the singularities" in such systems of (nonlinear) ODEs. Adaptive solvers typically rely on heuristic detection of increasing error estimates to decide on a switch in methods (or just decreasing step size). Feb 22 '15 at 19:35