# Set of integrators do not consistently solve an equation in Python

I must solve the following second order differential equation:

$$\delta \phi^{''}_{\mathbf{k}}+(3-\epsilon)\delta \phi^{'}_{\mathbf{k}}+\left(\frac{k^2}{a^2 H^2}+\frac{V_{,\phi\phi}}{H^2}-6\epsilon +4\epsilon\eta -2\epsilon^2\right)\delta \phi_{\mathbf{k}}=0;$$

which is that of scalar field perturbations at linear order in cosmic inflation. $$\epsilon$$, $$\eta$$, $$\phi$$, $$V_{,\phi\phi}$$, $$H$$ and $$a$$ denote background quantities which can be solved separately (i.e., their corresponding field equations do not contain field perturbation terms). $$k$$ represents the Fourier mode to which the field perturbation is ascribed and can be viewed as constant.

This equation can also be written as follows:

$$\frac{1}{a^3 \phi^{'2}H}\left[a^3 \phi^{'2} H \left(\frac{\delta \phi_{\mathbf{k}}}{\phi^{'}}\right)^{'}\right]^{'} =-\frac{k^2}{a^2H^2} \frac{\delta \phi_{\mathbf{k}}}{\phi^{'}}.$$

When I proceed to solve the first equation above, I use $$\delta \phi_{\mathbf{k}}$$ to solve the second one at the same time. However, it turns out that I do not get to solve the equations consistently (i.e., $$a^3 \phi^{'2} H \left(\frac{\delta \phi_{\mathbf{k}}}{\phi^{'}}\right)^{'}$$ obtained from the first equation does not agree with the one got from the last one, by far).

To solve the equations, I am employing various SciPy integrator solvers (RK45, RK23, LSODA, BDF, etc...), and none of them sorts out the discrepancy successfully.

Is there any suggestion about how I could get a consistent (more accurate) solution from those equations? Let me know if further details are needed.

• If you have two ways of solving the same equation, and they do not converge to the same answer as you make the time step small (or tolerance small), then at least one of them is wrongly implemented. Since once can assume that the SciPy integrators are correctly implemented, one might assume that you have a bug in your implementation. – Wolfgang Bangerth Aug 11 at 19:35
• In case the SciPy integrators were correctly implemented (and so were the equations in my code) but not accurate enough for my purposes, would you know of any other alternatives to solve these equations and get more consistent solutions? – J.J Aug 18 at 8:48
• I don't think the integrators are the problem. I'd suggest there's a bug in your implementation -- because 95% of the time, the problem is in user code. – Wolfgang Bangerth Aug 18 at 16:33