Skip to content

Adams-Moulton solvers are a family of implicit solvers that perform a single right-hand-side solve per time step and achieve high-order accuracy by storing the derivative at several previous time steps.

Adams-Moulton methods are multi-step solvers and thus require not only an initial condition, but also the solution at the first few time points depending on the order of the method. For example, AM3 requires the solution at times \(t_0, t_0+k, t_0+2k\) where \(t_0\) is the initial time and \(k\) is the temporal step-size.

This implementation assumes that you only have the solution at \(t_0\) and accepts another time-integrator as a seed. For example, we may use RK4 (a single step method). You must also decide the step-size of the single-step method such that the step-size of the Adams-Moulton solver is an integer multiple.

For example, 

    solver = AM4(RK4(), seed_steps_per_step=2)
Will use the initial condition and RK4 to take 8 time times steps each at half the step size. This will generate solutions at times \(t_0, t_0 + k/2, t_0 + k, t0 + 3k/2, ..., t_0 + 4k\). It will then sub-sample these as necessary to use as appropriate seed steps.

It is important to ensure that the order of the seed method and the number of seed steps per step are at least as accurate as the order of the Adams-Moulton method. For example 

solver = AM4(RK4(), seed_steps_per_step=1)
would not give a fifth order method because RK4 will generate seed steps up to fourth order accuracy. In contrast 
solver = AM5(RK4(), seed_steps_per_step=2)
will generate seed steps to eighth order accuracy, relative to the step size of the Adams-Moulton solver.

adams_moulton

A module containing Adams-Moulton solvers of several orders.

AdamsMoultonAbstract

Bases: TimeIntegrator

An abstract class for Adams-Moulton (AM) solvers. All AM solvers will inherit these methods. Additionally, they inherit the methods of time_integrator.TimeIntegrator as well.

Source code in src/odeiter/adams_moulton.py 
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
class AdamsMoultonAbstract(TimeIntegrator):
    """
    An abstract class for Adams-Moulton (AM) solvers.
    All AM solvers will inherit these methods. Additionally,
    they inherit the methods of [`time_integrator.TimeIntegrator`](time_integrator.md)
    as well.
    """

    def __init__(
        self,
        seed: TimeIntegrator,
        seed_steps_per_step: int,
        root_finder: RootFinder = DefaultRootFinder,
    ):
        """
        Parameters:
            seed: another time integrator used to take the first few steps.
            seed_steps_per_step: the number of seed steps taking per step
                of the AB integrator.
            root_finder: a function of the form `root_finder(f, u0)` that returns
                the solution to `f(u) = 0` using `u0` as an approximate soluiton.
        """
        self.seed = seed
        self.seed_steps_per_step = seed_steps_per_step
        self.seed_steps = self.order - 2
        self.root_finder = root_finder.solve

    @abstractproperty
    def name(self) -> str:
        """
        Returns:
            The name of the method
        """
        ...

    @abstractproperty
    def order(self) -> int:
        """
        Returns:
            The order of the method
        """
        ...

    @abstractproperty
    def fs_coeffs(self) -> list[float]:
        """
        Returns:
            The interpolation coefficients of the method.
        """
        ...

    def update(
        self,
        t: float,
        u: np.ndarray[float],
        rhs: Callable[[float, np.ndarray[float]], np.ndarray[float]],
        fs: np.ndarray[float],
        delta_t: float,
    ) -> np.ndarray[float]:
        """
        Compute the next time step. You probably want `solution_generator` instead.

        Parameters:
            t: The current time.
            u: The solution at the current time-step.
            rhs: The right-hand-side of the system as a function `rhs(t, u) -> u'`.
            fs: the right-hand-side at several previous time-steps.
            delta_t: the temporal step-size.

        Returns:
            The solution at the next time step.
        """
        const_vec = u + delta_t * sum(
            c * f for c, f in zip(self.fs_coeffs[-1:0:-1], fs)
        )

        def func(x):
            return -x + const_vec + delta_t * self.fs_coeffs[0] * rhs(t + delta_t, x)

        return self.root_finder(func, u)

    def solution_generator(
        self,
        u0: np.ndarray[float],
        rhs: Callable[[float, np.ndarray[float]], np.ndarray[float]],
        time: TimeDomain,
    ) -> Generator[np.ndarray[float], None, None]:
        """Create a generator that yields the solution for each time in `time`.

        Parameters:
            u0: The initial condition of the system.
                Must be the same size as the system.
            rhs: The right-hand-side as a function with signature `rhs(t, u) -> u'`.
            time: The discretized time domain from.

        Returns:
            A generator that yields the solution at each time in `time.array`.
        """
        seed_steps = self.order - 1
        seed_time = TimeRay(
            time.start,
            time.spacing / self.seed_steps_per_step,
        )
        t = time.start
        fs = deque(maxlen=self.order - 1)
        for u in islice(
            self.seed.solution_generator(u0, rhs, seed_time),
            0,
            self.seed_steps_per_step * seed_steps + 1,
            self.seed_steps_per_step,
        ):
            yield u
            fs.append(rhs(t, u))
            t += time.spacing

        for t in time.array[seed_steps:-1]:
            u = self.update(t, u, rhs, fs, time.spacing)
            yield u
            fs.append(rhs(t + time.spacing, u))

__init__(seed, seed_steps_per_step, root_finder=DefaultRootFinder)

Parameters:

Name Type Description Default
seed TimeIntegrator

another time integrator used to take the first few steps.

 required
seed_steps_per_step int

the number of seed steps taking per step of the AB integrator.

 required
root_finder RootFinder

a function of the form root_finder(f, u0) that returns the solution to f(u) = 0 using u0 as an approximate soluiton.

 DefaultRootFinder
Source code in src/odeiter/adams_moulton.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def __init__(
    self,
    seed: TimeIntegrator,
    seed_steps_per_step: int,
    root_finder: RootFinder = DefaultRootFinder,
):
    """
    Parameters:
        seed: another time integrator used to take the first few steps.
        seed_steps_per_step: the number of seed steps taking per step
            of the AB integrator.
        root_finder: a function of the form `root_finder(f, u0)` that returns
            the solution to `f(u) = 0` using `u0` as an approximate soluiton.
    """
    self.seed = seed
    self.seed_steps_per_step = seed_steps_per_step
    self.seed_steps = self.order - 2
    self.root_finder = root_finder.solve

name()

Returns:

Type Description
str

The name of the method
Source code in src/odeiter/adams_moulton.py 
40
41
42
43
44
45
46
@abstractproperty
def name(self) -> str:
    """
    Returns:
        The name of the method
    """
    ...

order()

Returns:

Type Description
int

The order of the method
Source code in src/odeiter/adams_moulton.py 
48
49
50
51
52
53
54
@abstractproperty
def order(self) -> int:
    """
    Returns:
        The order of the method
    """
    ...

fs_coeffs()

Returns:

Type Description
list[float]

The interpolation coefficients of the method.
Source code in src/odeiter/adams_moulton.py 
56
57
58
59
60
61
62
@abstractproperty
def fs_coeffs(self) -> list[float]:
    """
    Returns:
        The interpolation coefficients of the method.
    """
    ...

update(t, u, rhs, fs, delta_t)

Compute the next time step. You probably want solution_generator instead.

Parameters:

Name Type Description Default
t float

The current time.

 required
u ndarray[float]

The solution at the current time-step.

 required
rhs Callable[[float, ndarray[float]], ndarray[float]]

The right-hand-side of the system as a function rhs(t, u) -> u'.

 required
fs ndarray[float]

the right-hand-side at several previous time-steps.

 required
delta_t float

the temporal step-size.

 required

Returns:

Type Description
ndarray[float]

The solution at the next time step.
Source code in src/odeiter/adams_moulton.py 
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def update(
    self,
    t: float,
    u: np.ndarray[float],
    rhs: Callable[[float, np.ndarray[float]], np.ndarray[float]],
    fs: np.ndarray[float],
    delta_t: float,
) -> np.ndarray[float]:
    """
    Compute the next time step. You probably want `solution_generator` instead.

    Parameters:
        t: The current time.
        u: The solution at the current time-step.
        rhs: The right-hand-side of the system as a function `rhs(t, u) -> u'`.
        fs: the right-hand-side at several previous time-steps.
        delta_t: the temporal step-size.

    Returns:
        The solution at the next time step.
    """
    const_vec = u + delta_t * sum(
        c * f for c, f in zip(self.fs_coeffs[-1:0:-1], fs)
    )

    def func(x):
        return -x + const_vec + delta_t * self.fs_coeffs[0] * rhs(t + delta_t, x)

    return self.root_finder(func, u)

solution_generator(u0, rhs, time)

Create a generator that yields the solution for each time in time.

Parameters:

Name Type Description Default
u0 ndarray[float]

The initial condition of the system. Must be the same size as the system.

 required
rhs Callable[[float, ndarray[float]], ndarray[float]]

The right-hand-side as a function with signature rhs(t, u) -> u'.

 required
time TimeDomain

The discretized time domain from.

 required

Returns:

Type Description
None

A generator that yields the solution at each time in time.array.
Source code in src/odeiter/adams_moulton.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def solution_generator(
    self,
    u0: np.ndarray[float],
    rhs: Callable[[float, np.ndarray[float]], np.ndarray[float]],
    time: TimeDomain,
) -> Generator[np.ndarray[float], None, None]:
    """Create a generator that yields the solution for each time in `time`.

    Parameters:
        u0: The initial condition of the system.
            Must be the same size as the system.
        rhs: The right-hand-side as a function with signature `rhs(t, u) -> u'`.
        time: The discretized time domain from.

    Returns:
        A generator that yields the solution at each time in `time.array`.
    """
    seed_steps = self.order - 1
    seed_time = TimeRay(
        time.start,
        time.spacing / self.seed_steps_per_step,
    )
    t = time.start
    fs = deque(maxlen=self.order - 1)
    for u in islice(
        self.seed.solution_generator(u0, rhs, seed_time),
        0,
        self.seed_steps_per_step * seed_steps + 1,
        self.seed_steps_per_step,
    ):
        yield u
        fs.append(rhs(t, u))
        t += time.spacing

    for t in time.array[seed_steps:-1]:
        u = self.update(t, u, rhs, fs, time.spacing)
        yield u
        fs.append(rhs(t + time.spacing, u))

AM2

Bases: AdamsMoultonAbstract

Adams Moulton of order 3.

Source code in src/odeiter/adams_moulton.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
class AM2(AdamsMoultonAbstract):
    """Adams Moulton of order 3."""

    @property
    def order(self):
        return 3

    @property
    def name(self):
        return f"AM2 (seed: {self.seed.name})"

    @property
    def fs_coeffs(self) -> list[float]:
        return [5 / 12, 8 / 12, -1 / 12]

AM3

Bases: AdamsMoultonAbstract

Adams Moulton of order 4.

Source code in src/odeiter/adams_moulton.py 
153
154
155
156
157
158
159
160
161
162
163
164
165
166
class AM3(AdamsMoultonAbstract):
    """Adams Moulton of order 4."""

    @property
    def order(self):
        return 4

    @property
    def name(self):
        return f"AM3 (seed: {self.seed.name})"

    @property
    def fs_coeffs(self) -> list[float]:
        return [9 / 24, 19 / 24, -5 / 24, 1 / 24]

AM4

Bases: AdamsMoultonAbstract

Adams Moulton of order 5.

Source code in src/odeiter/adams_moulton.py 
169
170
171
172
173
174
175
176
177
178
179
180
181
182
class AM4(AdamsMoultonAbstract):
    """Adams Moulton of order 5."""

    @property
    def order(self):
        return 5

    @property
    def name(self):
        return f"AM4 (seed: {self.seed.name})"

    @property
    def fs_coeffs(self) -> list[float]:
        return [251 / 720, 646 / 720, -264 / 720, 106 / 720, -19 / 720]