I'm trying to numerical solve a 2nd order differential equation, I know the usual way of doing this is using Runge-Kutta's Method. But then I wonder if I could do it with Adam-Bashford and Adam-Bashford predictor-corrector. To my surprise, although theoretically, it is possible to do so, I didn't end up finding any implementation nor straight forward formula, so I wrote this:
# AB
Y[i + 1] = Y[i] + h*(55.0 * U[i] - 59.0 * U[i-1] + 37.0 * U[i-2] - 9.0 * U[i-3])/24.0
U[i + 1] = U[i] + h*(55.0 * dudx(X[i],Y[i],U[i]) - 59.0 * dudx(X[i-1],Y[i-1],U[i-1]) + 37.0 * dudx(X[i-2],Y[i-2],U[i-2]) - 9.0 * dudx(X[i-3],Y[i-3],U[i-1]))/24.0
#AB predictor-corrector
Y[i + 1] = Y[i] + h*(55.0 * U[i] - 59.0 * U[i-1] + 37.0 * U[i-2] - 9.0 * U[i-3])/24.0
U[i + 1] = U[i] + h*(55.0 * dudx(X[i],Y[i],U[i]) - 59.0 * dudx(X[i-1],Y[i-1], U[i-1]) + 37.0 * dudx(X[i-2],Y[i-2], U[i-2]) - 9.0 * dudx(X[i-3],Y[i-3],U[i-3]))/24.0
Y[i + 1] = Y[i] + h*(9.0 * U[i+1] + 19.0 * U[i] - 5.0 * U[i-1] + U[i-2])/24.0
U[i + 1] = U[i] + h*(9.0 * dudx(X[i+1], Y[i + 1], U[i+1]) + 19.0 * dudx(X[i],Y[i], U[i]) - 5.0 * dudx(X[i-1],Y[i-1], U[i-1]) + dudx(X[i-2],Y[i-2], U[i-2]))/24.0
where $u=\frac{dy}{dx}$ and $dudx=\frac{d^2 y}{dx^2}$. I believed that this would result in a better estimate of the solution. But it seems AB is less accurate than rk4 and ABP-C is just a little better than rk4, Have I done something wrong? If not, Is there any way that someone can explain why this happens?
This is the result I got for the equation:
$$
\left\{\begin{array}{l}
y''+y'-6y=0\\
y(0)=1, y'(0)=2
\end{array}\right.
$$
where $y(x)=e^{2x}$ with this different methods:
Best Answer
There is one error and one non-standard feature in your formulas, and thus apparently also in your computations. In the AB4 formula for the
Ycomponent, the last index should bei-3as the ones directly before it, noti-1. In the repetition as predictor for AM4 you have that correct. In the corrector you are using the freshly computed value forY[i+1]in the computation ofU[i+1]. This will produce a difference to vector-based standard implementations. This is basically the difference between a Gauß-Seidel scheme and a Jacobi scheme. If one were to continue the corrector iteration until numerical convergence, there should be no differences above floating point noise.To check if you implemented the method correctly, you can compare the errors when dividing the step size by
N=2,4,8etc. You should observe a reduction in error by a factor $16$ in each step. For a more visual confirmation, divide the step size byN=10, 100, etc. Then the number of correct digits should increase by 4 in each step. I computed such tables below, they show that clearly in the columns forN=10. It should be true in general that the columns have about 2, 6 and 10 correct digits after the dot, respectively, RK4 and AM4 partially a little more, AB4 a little less.Runge-Kutta 4
Adams-Bashford 4
Adams-Moulton 4-1