Please correct if there a mistake on my Explicit formula of 7th Steps Adam-Bashforth-Moulton

Question

I'm searching in the internet, and i can't find about the formula (Sadly).

But if you know or found the amazing sites that gives some formula of Explicit Adam-Bashforth for any higer order (including 7th) in formal equation or code please share it to me.

And for my final choice, finally i found the formula by myself. I derive it manually with lagrange interpolation of order 6. (It's true isn't it? Because ABM 7th needs 7 points to constructs its formula).

Here is my formula.

Predictor of 7th ABM :

$$y^*(x_{n+1})=y_n+\frac{h}{60480}(198721f_n-894576 f_{n-1} +705549 f_{n-2}-688256 f_{n-3}+407139 f_{n-4}-134472 f_{n-5}+19087 f_{n-6})$$

Corrector of 7th ABM :

$$y(x_{n+1})=y_n+\frac{h}{60480}(65112 f_n-92922 f_{n-1} +37504 f_{n-2}-20211 f_{n-3}+6312 f_{n-4}-863 f_{n-5}-863 f_{n+1}^*)$$

Please verify my formula if there is a mistake. (On some numbers or anything).

Anyway, is it true when i'm using 6th Lagrange Interpolation with 7 points to construct the 7th Order (p=7) Adam-Bashforth-Moulton?

Thank in advance!

Answer 1

In the generic test case $y'=\lambda y$, the solution is $y=Ce^{λt}$. Inserting this exact solution into the AB scheme $$ y_{n+1}=y_n+h\sum_{k=0}^{p}b_kf(t_{n-k},y_{n-k}) $$ in order to compute the truncation error, and cancelling common factors results in the formula $$ e^{λh}=1+λh\sum_{k=0}^{p}b_ke^{-kλh}+O(h^{p+1}). $$ So the aim is to approximate $\frac{e^{λh}-1}{λh}$ as a polynomial in $e^{-λh}$ around $h=0$. Selecting another variable $z=1-e^{-λh}$ with $h=0\iff z=0$, the polynomial in $e^{-λh}$ transforms into a polynomial $q(z)$ in $z$ with likewise unknown coefficients. These are determined by the equation in truncated power series $$ z=-\log(1-z)(1-z)q(z)+O(z^{p+1}). $$ Good computer algebra systems enable working with truncated power series, for instance in CAS Magma (online calculator) one can execute this via

PS<z>:=PowerSeriesRing(Rationals());
for p in [2..9] do
    // solve the power series equation
    q:=-z/((1-z)*Log(1-z+O(z^(p+3))));
    // extract the polynomial as polynomial
    q:=Truncate(q+O(z^p));
    // shift so that the coefficients are the b_k
    c:=Coefficients(Evaluate(q, 1-z));
    // for cosmetic reasons use a common denominator in the fractions
    den := LCM([Denominator(cc): cc in c]);
    Sprintf("%o: %o/%o", p, [ den*cc: cc in c],den);
end for;

giving the table of right sides

2: [ 3, -1 ]/2
3: [ 23, -16, 5 ]/12
4: [ 55, -59, 37, -9 ]/24
5: [ 1901, -2774, 2616, -1274, 251 ]/720
6: [ 4277, -7923, 9982, -7298, 2877, -475 ]/1440
7: [ 198721, -447288, 705549, -688256, 407139, -134472, 19087 ]/60480
8: [ 434241, -1152169, 2183877, -2664477, 2102243, -1041723, 295767, -36799 ]/120960
9: [ 14097247, -43125206, 95476786, -139855262, 137968480, -91172642, 38833486, -9664106, 1070017 ]/3628800

For Adams-Moulton $$ y_{n+1}=y_n+h\sum_{k=0}^{p}b_kf(t_{n+1-k},y_{n+1-k}) $$ you get that the highest indices on both sides are the same, so the test case results, after cancelling one power more, in the equation $$ 1=e^{-λh}+λh\sum_{k=0}^{p}b_ke^{-kλh}+O(h^{p+1}). $$ That is, the task is now to approximate $\frac{1-e^{-λh}}{λh}$ with a polynomial in $e^{-λh}$ with remainder $O(h^p)$, that is, similarly to the above from

PS<z>:=PowerSeriesRing(Rationals());
for p in [2..9] do
    q:=-z/(Log(1-z+O(z^(p+3))));
    q:=Truncate(q+O(z^p));
    c:=Coefficients(Evaluate(q, 1-z));
    den := LCM([Denominator(cc): cc in c]);
    Sprintf("%o: %o/%o", p, [ den*cc: cc in c],den);
end for;

the table

2: [ 1, 1 ]/2
3: [ 5, 8, -1 ]/12
4: [ 9, 19, -5, 1 ]/24
5: [ 251, 646, -264, 106, -19 ]/720
6: [ 475, 1427, -798, 482, -173, 27 ]/1440
7: [ 19087, 65112, -46461, 37504, -20211, 6312, -863 ]/60480
8: [ 36799, 139849, -121797, 123133, -88547, 41499, -11351, 1375 ]/120960
9: [ 1070017, 4467094, -4604594, 5595358, -5033120, 3146338, -1291214, 312874, -33953 ]/3628800

Extracting specifically for order $p=7$, the coefficients are thus

AB: [ 198721, -447288, 705549, -688256, 407139, -134472, 19087 ]/60480
AM: [ 19087, 65112, -46461, 37504, -20211, 6312, -863 ]/60480

which show differences in several places to your result.