Skip to content

Backward Differentiation (BDF) 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 solution at several previous time steps.

BDF 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, BDF3 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 = BDF5(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 BDF method. For example 

solver = BDF5(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 = BDF5(RK4(), seed_steps_per_step=2)
will generate seed steps to eighth order accuracy, relative to the step size of the BDF solver.

backward_differentiation

BackwardDifferentiationAbstract

Bases: TimeIntegrator

An abstract class for Backward Differentiation (BDF) solvers. All BDF solvers will inherit these methods. Additionally, they inherit the methods of time_integrator.TimeIntegrator as well.

Source code in src/odeiter/backward_differentiation.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
132
133
134
class BackwardDifferentiationAbstract(TimeIntegrator):
    """
    An abstract class for Backward Differentiation (BDF) solvers.
    All BDF 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 - 1
        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 us_coeffs(self) -> list[float]:
        """
        Returns:
            The solution interpolation coefficients of the method.
        """
        ...

    @abstractproperty
    def f_coeff(self) -> float:
        """
        Returns:
            The derivative interpolation coefficient of the method.
        """
        ...

    def update(
        self,
        t: float,
        us: np.ndarray[float],
        rhs: Callable[[float, np.ndarray[float]], 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.
            us: The solution at the current time-step and several previous time-steps.
            rhs: The right-hand-side of the system as a function `rhs(t, u) -> u'`.
            delta_t: the temporal step-size.

        Returns:
            The solution at the next time step.
        """
        const_vec = sum(c * u for c, u in zip(self.us_coeffs[::-1], us))

        def func(x):
            return -x + const_vec + delta_t * self.f_coeff * rhs(t + delta_t, x)

        return self.root_finder(func, us[-1])

    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
        us = deque(maxlen=self.order)
        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
            us.append(u)

        for t in time.array[seed_steps:-1]:
            u = self.update(t, us, rhs, time.spacing)
            yield u
            us.append(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/backward_differentiation.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 - 1
    self.root_finder = root_finder.solve

name()

Returns:

Type Description
str

The name of the method
Source code in src/odeiter/backward_differentiation.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/backward_differentiation.py 
48
49
50
51
52
53
54
@abstractproperty
def order(self) -> int:
    """
    Returns:
        The order of the method
    """
    ...

us_coeffs()

Returns:

Type Description
list[float]

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

f_coeff()

Returns:

Type Description
float

The derivative interpolation coefficient of the method.
Source code in src/odeiter/backward_differentiation.py 
64
65
66
67
68
69
70
@abstractproperty
def f_coeff(self) -> float:
    """
    Returns:
        The derivative interpolation coefficient of the method.
    """
    ...

update(t, us, rhs, delta_t)

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

Parameters:

Name Type Description Default
t float

The current time.

 required
us ndarray[float]

The solution at the current time-step and several previous time-steps.

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

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

 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/backward_differentiation.py 
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
def update(
    self,
    t: float,
    us: np.ndarray[float],
    rhs: Callable[[float, np.ndarray[float]], 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.
        us: The solution at the current time-step and several previous time-steps.
        rhs: The right-hand-side of the system as a function `rhs(t, u) -> u'`.
        delta_t: the temporal step-size.

    Returns:
        The solution at the next time step.
    """
    const_vec = sum(c * u for c, u in zip(self.us_coeffs[::-1], us))

    def func(x):
        return -x + const_vec + delta_t * self.f_coeff * rhs(t + delta_t, x)

    return self.root_finder(func, us[-1])

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/backward_differentiation.py 
 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
132
133
134
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
    us = deque(maxlen=self.order)
    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
        us.append(u)

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

BDF2

Bases: BackwardDifferentiationAbstract

Backward Differentiation of order 2.

Source code in src/odeiter/backward_differentiation.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
class BDF2(BackwardDifferentiationAbstract):
    """Backward Differentiation of order 2."""

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

    @property
    def order(self):
        return 2

    @property
    def us_coeffs(self) -> list[float]:
        return [4 / 3, -1 / 3]

    @property
    def f_coeff(self) -> float:
        return 2 / 3

BDF3

Bases: BackwardDifferentiationAbstract

Backward Differentiation of order 3.

Source code in src/odeiter/backward_differentiation.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
class BDF3(BackwardDifferentiationAbstract):
    """Backward Differentiation of order 3."""

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

    @property
    def order(self):
        return 3

    @property
    def us_coeffs(self) -> list[float]:
        return [18 / 11, -9 / 11, 2 / 11]

    @property
    def f_coeff(self) -> float:
        return 6 / 11

BDF4

Bases: BackwardDifferentiationAbstract

Backward Differentiation of order 4.

Source code in src/odeiter/backward_differentiation.py 
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
class BDF4(BackwardDifferentiationAbstract):
    """Backward Differentiation of order 4."""

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

    @property
    def order(self):
        return 4

    @property
    def us_coeffs(self) -> list[float]:
        return [48 / 25, -36 / 25, 16 / 25, -3 / 25]

    @property
    def f_coeff(self) -> float:
        return 12 / 25

BDF5

Bases: BackwardDifferentiationAbstract

Backward Differentiation of order 5.

Source code in src/odeiter/backward_differentiation.py 
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
class BDF5(BackwardDifferentiationAbstract):
    """Backward Differentiation of order 5."""

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

    @property
    def order(self):
        return 5

    @property
    def us_coeffs(self) -> list[float]:
        return [300 / 137, -300 / 137, 200 / 137, -75 / 137, 12 / 137]

    @property
    def f_coeff(self) -> float:
        return 60 / 137