In this tutorial you’ll learn how to

  • create custom integration schemes

Custom Integration Schemes

In this tutorial we want to estimate \(\pi\) with the following equation:

\(\pi = 4 \int_0^1 \sqrt{1-t^2}\mathrm{d}t\)

We set up the simulation frame as explained in the previous tutorials.

[1]:
from simframe import Frame
[2]:
sim = Frame()

Adding field for :math:`pi` and integration variable

[3]:
sim.addfield("pi", 0., description="Approximation of pi")
sim.addintegrationvariable("t", 0.)

Differentiator

[4]:
import numpy as np

def f(frame, x, Y):
    return 4.*np.sqrt(1-x**2)
[5]:
sim.pi.differentiator = f

Step size

We set the step size to 0.25, i.e., the integral function is only evaluated four times in the simulation.

[6]:
def dt(frame):
    return 0.25
[7]:
sim.t.updater = dt

We do not need to write outputs for this model. We only have to tell the integrator when to stop the simulation, i.e., the upper bound of the integral.

[8]:
sim.t.snapshots = [1.]

Setting the integrator

[9]:
from simframe import Integrator
from simframe import Instruction
from simframe import schemes
[10]:
sim.integrator = Integrator(sim.t)
[11]:
sim.integrator.instructions = [Instruction(schemes.expl_1_euler, sim.pi)]

Running the simulation

[12]:
sim.run()
Execution time: 0:00:00

Results

[13]:
from IPython.display import Markdown as md
def print_table(sim):
    return md("| |$\pi$|rel. error|\n|-|-|-|\n|real|{:10.8f}||\n|approx.|{:10.8f}|{:9.3e}|".format(np.pi,sim.pi,np.abs(np.pi-sim.pi)/np.pi))
[14]:
print_table(sim)
[14]:

\(\pi\)

rel. error

real

3.14159265

approx.

3.49570907

1.127e-01

The relative error is about 11 %.

There are in principle two way how we can improve this estimate. One is to decrease the stepsize. You can rerun the notebook with a smaller step size and compare the results.

Another method is to use more sophisticated integration schemes. We are now going to implement the 2nd-order Midpoint Method and the 4th-order Runge-Kutta method.

Writing custom integration schemes

All we have to do is write a function that takes the current value of the integration variable x0, the current value of the variable to be integrated Y0, and the step size dx. The function needs to return the change dY of our variable Y0 after the integration. We can access the derivative of the variable with Y0.derivative(x, Y).

For the midpoint method this looks as follows.

[15]:
def midpoint(x0, Y0, dx, *args, **kwargs):
    x1 = x0 + 0.5*dx
    Y1 = Y0 + 0.5*dx*Y0.derivative(x0, Y0)
    return dx*Y0.derivative(x1, Y1)

All we have to do now is to create an integration scheme from this function. This can be done by using AbstractScheme provided by simframe.

[16]:
from simframe.integration import Scheme
[17]:
expl_2_midpoint = Scheme(midpoint)

expl_2_midpoint is now our new integration scheme. The naming convention is <expl/impl>_<order>_<name><_other>.

We can now assign a new instruction set using our new method to the integrator just as with the 1st-order Euler method.

[18]:
sim.integrator.instructions = [Instruction(expl_2_midpoint, sim.pi)]

Before we restart the simulation we have to reset to the initial conditions.

[19]:
sim.t = 0
sim.pi = 0
[20]:
sim.run()
Execution time: 0:00:00
[21]:
print_table(sim)
[21]:

\(\pi\)

rel. error

real

3.14159265

approx.

3.18392922

1.348e-02

The error is now reduced to 1.3 % only by using a higher order method.

Note: the higher order method needs more operations and is therefore slower. This does not really matter in this case, but might be important for more complex simulations.

4th-order Runge-Kutta

The scheme function for the 4th-order explicit Runge-Kutta method looks as follows.

[22]:
def rk4(x0, Y0, dx, *args, **kwargs):
    k1 = Y0.derivative(x0         , Y0            )
    k2 = Y0.derivative(x0 + 0.5*dx, Y0 + 0.5*dx*k1)
    k3 = Y0.derivative(x0 + 0.5*dx, Y0 + 0.5*dx*k2)
    k4 = Y0.derivative(x0 +     dx, Y0 +     dx*k3)
    return dx*(1/6*k1 + 1/3*k2 + 1/3*k3 + 1/6*k4)
[23]:
expl_4_rungekutta = Scheme(rk4)
[24]:
sim.integrator.instructions = [Instruction(rk4, sim.pi)]
[25]:
sim.t = 0
sim.pi = 0
[26]:
sim.run()
Execution time: 0:00:00
[27]:
print_table(sim)
[27]:

\(\pi\)

rel. error

real

3.14159265

approx.

3.12118917

6.495e-03

With this method the error is reduced down to 0.6 %.

Available integration schemes

But before you take effort into developing your own integration schemes, take a look at the schemes already provided by simframe:

[28]:
_  = "|Scheme|Description|\n"
_ += "|------|-----------|\n"
for s in schemes.__dir__():
    if s == "update": continue
    if not s.startswith("_"): _ += "|"+s+"|"+schemes.__dict__[s].description+"|\n"
md(_)
[28]:

Scheme

Description

expl_1_euler

Explicit 1st-order Euler method

expl_2_fehlberg_adptv

Explicit adaptive 2nd-order Fehlberg’s method

expl_2_heun

Explicit 2nd-order Heun’s method

expl_2_heun_euler_adptv

Explicit adaptive 2nd-order Heun-Euler method

expl_2_midpoint

Explicit 2nd-order midpoint method

expl_2_ralston

Explicit 2nd-order Ralston’s method

expl_3_bogacki_shampine_adptv

Explicit adaptive 3rd-order Bogacki-Shampine method

expl_3_gottlieb_shu_adptv

Explicit adaptive 3rd-order Gottlieb-Shu method

expl_3_heun

Explicit 3rd-order Heun’s method

expl_3_kutta

Explicit 3rd-order Kutta’s method

expl_3_ralston

Explicit 3rd-order Ralston’s method

expl_3_ssprk

Explicit 3rd-order Strong Stability Preserving Runge-Kutta method

expl_4_38rule

Explicit 4th-order 3/8 rule method

expl_4_ralston

Explicit 4th-order Ralston’s method

expl_4_runge_kutta

Explicit 4th-order classical Runge-Kutta method

expl_5_cash_karp_adptv

Explicit adaptive 5th-order Cash-Karp method

impl_1_euler_direct

Implicit 1st-order direct Euler method

impl_1_euler_gmres

Implicit 1st-order Euler method with GMRES solver

impl_2_midpoint_direct

Implicit 2nd-order direct midpoint method