4. Custom Integration Schemes¶
In this tutorial we want to estimate \(\pi\) with the following equation:
\(\pi = 4 \int\limits_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 adaptive 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 |
expl_5_dormand_prince_adptv |
Explicit adaptive 5th-order Dormand-Prince 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 |