# Calculating the bounce solution numerically

I would like to obtain a numerical solution to the following example bounce equations,

\begin{align*} \frac{\partial^2 a}{\partial t^2} &= \frac{1}{t^2} a(1-a)(1-3a)-\frac{b^2}{2}(1-a)\\ \frac{\partial^2 b}{\partial t^2} &= \frac{1}{t^2}b(1-a)^2-\frac{1}{t}\frac{\partial b}{\partial t}+b(b^2-1) \end{align*}

with formal boundary conditions

$$a(0)=b(0)=0, \; a(\infty)=b(\infty)=1,$$

although $$a(x)=b(x)=1$$ for any large value of $$x$$ should be fine as well. If possible, I would like to do so by largely using built-in Mathematica functions. I have changed some of the numerical prefactors, but the essentials of the system are the same. Whenever I carry out this calculation with NDSolve, Mathematica complains about the system being stiff - I suspect this is due to the $$\frac{1}{t^{(2)}}$$ factors evaluated close to the origin. As a beginner to Mathematica and numerical methods, I really have no further tricks up my sleeve at the moment, so I would appreciate any good hint. I should mention that I have been able to find a solution by implementing a particular numerical method (an elaborate version of the Newton algorithm), but it would be much easier, arguably less error-prone and hopefully better scalable to more equations if I found a way of having Mathematica take care of it largely by itself. Thank you very much.

The problem you are running into is that $$t = 0$$ is a singular point of the system. This is (as you suspect) due to the factors of $$1/t$$ and $$1/t^2$$ in the equations. A straightforward implementation of an ODE solver will not be able to solve all the way to $$t = 0$$ in general, since the highest-order derivatives will diverge at that point for a general solution.

One way to get around this for a singular system is to use Mathematica to find an approximate solution near the singular point. This allows us to "move the boundary" away from the singular point and solve over a domain $$[\epsilon, T]$$ (for some large $$T$$) that doesn't contain any points where Mathematica will obviously have trouble. By finding an approximate solution near $$t = 0$$, we can find the approximate value of the solution at $$t = \epsilon > 0$$, perhaps in terms of some arbitrary parameters. If $$\epsilon$$ is sufficiently close to zero, then $$a(\epsilon)$$ & $$b(\epsilon)$$ will be "close" to their "true" values in a solution with $$a(0) = b(0) = 0$$, and solving the equations over the domain $$t \in [\epsilon, T]$$ (for some large $$T$$) will give you a solution that should be "close" to the true solution.

In practical terms: The solution should (hopefully) be expressible as a power series in terms of $$t$$; and so we try an ansatz of a series solution that goes to zero at there.

eqns = {a''[t] - (1/t^2 a[t] (1 - a[t]) (1 - 3 a[t]) - b[t]^2/2 (1 - a[t])),
b''[t] - (1/t^2 b[t] (1 - a[t])^2 - 1/t b'[t] + b[t] (b[t]^2 - 1))}
seriesrules = {a -> Function[t, SeriesData[t, 0, {a1, a2, a3, a4, a5}, 1, 6, 1]],
b -> Function[t, SeriesData[t, 0, {b1, b2, b3, b4, b5}, 1, 6, 1]]}
coeffrules = SolveAlways[Normal[eqns /. seriesrules] == 0, t]

(* {{a5 -> 0, b5 -> (b1 (11 + 96 b1^2))/2112, b4 -> 0, a4 -> -(b1^2/22),
b3 -> -(b1/8), a3 -> 0, b2 -> 0, a2 -> 0, a1 -> 0}} *)


This implies that near $$r = 0$$, the solution will have the form

approxsoln = Normal[({a[eps], b[eps]} /. seriesrules) /. coeffrules]

(* {-(1/22) b1^2 eps^4, b1 eps - (b1 eps^3)/8 + (b1 (11 + 96 b1^2) eps^5)/2112} *)


One can in principle then put the boundary at eps and integrate out to some large T. By tweaking the value of b1, one can quickly narrow down the value of b1 that leads to boundary condition one wants. (This may be equivalent to the "elaborate Newton algorithm" method you describe in your question.)

Alternately, one could use b1 and eps as parameters in a ParametricNDSolve call, and then use the root-finding methods of FindRoot to find the value for which {a[inf], b[inf]} takes on a particular value.

This said, I suspect that for the particular equations you've provided, there is no solution that is analytic at $$t = 0$$, for the simple solution that there is only one free parameter in the series expansion. This means you effectively have to "hit a target" in 2-D ($$a(\infty) = b(\infty) = 1$$) with only one degree of freedom; it is likely that there is no value of $$b_1$$ that leads to such a solution. Perhaps, though, the techniques I've described here will be useful for your actual problem.

• A couple of things. First of all, thank you for your answer. Secondly, I do have a couple of questions on what you suggested. What exactly do you mean by ''putting the boundary at eps''? I understand that if I set eps=0, that series expansion of the functions will fulfill the BC at the origin, if that is what you mean; then I'd be left with setting a[x]=b[x]=1 for some arbitrarily large x in order to obtain b1. Or do you rather mean I set the left boundary at a symbolic eps and integrate the original equations up to some large t? I appreciate your effort, I'm afraid I'm just a bit lost. – Al Waurora Jun 4 '20 at 22:20
• Also, to get it out of the way - my original solution involves a variant of the finite difference method, which in the end leaves me with a discretized system of algebraic equations of which the Newton algorithm takes care. My algorithm works fine for this system of equations, but it does not yield the results I expect after adding another function. Today I've discovered that v12 of Mathematica has implemented in NDSolve the Method->FiniteElement (which I understand is their implementation of the Finite Difference method), and I was hoping someone could perhaps shed some light onto it. – Al Waurora Jun 4 '20 at 22:28
• @AlWaurora: I've edited the first couple of paragraphs to better explain the technique. – Michael Seifert Jun 4 '20 at 22:40
• Again, thank you very much. – Al Waurora Jun 5 '20 at 8:29