Scheme

The Scheme class may be used to easily create new schemes.

Let us examine the code to create the Explicit Euler scheme:

from odelab import Scheme

class ExplicitEuler(Scheme):
    def delta(self, t, u0, h):
        return h, h*self.system.f(t,u0)

The method delta() takes as argument the current time t and state u0, as well as the time step h. It returns a time increment and a state increment.

An object of the class Scheme has access to the attribute system, which gives access to the system at hand.

If the method is implicit, all you have to do is to implement the get_residual() method instead. Here is an example showing how to implement the Implicit Euler scheme:

class ImplicitEuler(Scheme):
    def get_residual(self, t, u0, h):
        def residual(du):
            return du - h*self.system.f(t+h, u0+du)
        return residual

Main Methods

The odelab.scheme.Scheme class provides several methods to assist the creation of numerical schemes.

Explicit Schemes

The method delta() may be used to directly provide the difference \(u_1 - u_0\). This is typically useful for explicit scheme. See for instance odelab.scheme.classic.RungeKutta4, or odelab.scheme.classic.ExplicitEuler.

Implicit Schemes

General Residual

For implicit schemes, the method get_residual() simply returns the residual of a function which root is the difference \(u_1 - u_0\). In other words, the residual is a function \(F\) such that

\[F(u_1 - u_0) = 0\]

Custom Residual

Sometimes, the residual’s root is not the difference \(u_1 - u_0\), and odelab supports that as well. The best way to understand how it works is to look at the implementation of delta(), which uses get_residual().

In short, the method delta runs a nonlinear solver on the residual provided by get_residual, using a guess provided by get_guess(). The root is then passed to the method reconstruct(), which job is to construct the difference \(u_1 - u_0\) from the root of the residual.

This flexibility is useful, for instance, if only one part of the method is implicit.

Time Adaptivity

The methods delta(), get_residual() or step() have h as a parameter, but it is just a convenience. The time step parameter is in fact stored by the Scheme object itself. A time adaptive scheme may thus change the next time step, or even implement advanced time-step control.

An example of time-adaptive scheme is given by RungeKutta34.

Initialization and External Solvers

It is the user’s responsibility to initialize the scheme object appropriately. Usually, the only necessary argument is the time-step:

ie = ImplicitEuler(h=.1)

Note that the solver calls the initialize() method before running the simulation. At this stage, the solver object has been initialized, and in particular has a initial condition and a System object. This is especially useful for using external numerical schemes, as those often require the right hand side and initial conditions to be known at initialization. An example of scheme using an external numerical scheme is given by ode15s, which uses SciPy’s function odeint.

Family of Numerical Schemes

odelab comes with a variety of numerical schemes to be used out of the box. For instance:

• generic Runge-Kutta methods
• generic linear multi-step methods
• generic explicit general linear methods
• generic exponential methods

”Classical“ Numerical Scheme

A collection of classical numerical scheme described in all the basic textbooks on differential equations.

class odelab.scheme.classic.ExplicitEuler(h=None, *args, **kwargs)

Explicit version of the Euler method, defined by:

\[u_1 = u_1 + hf(t_0, u_0)\]

class odelab.scheme.classic.ImplicitEuler(h=None, *args, **kwargs)

The standard implicit Euler scheme, defined by:

\[u_1 = u_0 + hf(t_1, u_1)\]

class odelab.scheme.classic.RungeKutta34(h=None, *args, **kwargs)

Adaptive Runge-Kutta of order four.

class odelab.scheme.classic.RungeKutta4(h=None, *args, **kwargs)

Standard Runge-Kutta of order 4.

The Scheme class

class odelab.scheme.Scheme(h=None, *args, **kwargs)

General Scheme class. Subclass this class to define a specific integration method.

adjust_stepsize(error)

Change the step size based on error estimation. To be overridden for a variable step size method.

delta(t, u0, h)

Compute the difference between current and next state.

get_guess(t, u0, h)

Default guess for the Newton iterations, assuming that the residual has the same size as u.

get_residual(t, u0, h)

Return the residual, which root is u1 - u0.

initialize(events)

Called the first time the scheme is used during a simulation.

reconstruct(root, t, u0, h)

Default reconstruction function. It assumes that the root is already delta u.

root_solver

alias of FSolve

step(t, u0, h)

Implementation of the Compensated Summation algorithm as described in [HaLuWa06] §VIII.5.

[HaLuWa06]

Hairer, Lubich, Wanner Geometric Numerical Integration, Springer, 2006.

class odelab.scheme.ode15s(**kwargs)

Simulation of matlab’s ode15s solver. It is a BDF method of max order 5