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{1t^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(1x**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\nreal{:10.8f}\napprox.{:10.8f}{:9.3e}".format(np.pi,sim.pi,np.abs(np.pisim.pi)/np.pi))
[14]:
print_table(sim)
[14]:
\(\pi\) 
rel. error 


real 
3.14159265 

approx. 
3.49570907 
1.127e01 
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 2ndorder Midpoint Method and the 4thorder RungeKutta 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 1storder 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.348e02 
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.
4thorder RungeKutta
The scheme function for the 4thorder explicit RungeKutta 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.495e03 
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]:
_ = "SchemeDescription\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 1storder Euler method 
expl_2_fehlberg_adptv 
Explicit adaptive 2ndorder Fehlberg’s method 
expl_2_heun 
Explicit 2ndorder Heun’s method 
expl_2_heun_euler_adptv 
Explicit adaptive 2ndorder HeunEuler method 
expl_2_midpoint 
Explicit 2ndorder midpoint method 
expl_2_ralston 
Explicit 2ndorder Ralston’s method 
expl_3_bogacki_shampine_adptv 
Explicit adaptive 3rdorder BogackiShampine method 
expl_3_gottlieb_shu_adptv 
Explicit adaptive 3rdorder GottliebShu method 
expl_3_heun 
Explicit 3rdorder Heun’s method 
expl_3_kutta 
Explicit 3rdorder Kutta’s method 
expl_3_ralston 
Explicit 3rdorder Ralston’s method 
expl_3_ssprk 
Explicit 3rdorder Strong Stability Preserving RungeKutta method 
expl_4_38rule 
Explicit 4thorder 3/8 rule method 
expl_4_ralston 
Explicit 4thorder Ralston’s method 
expl_4_runge_kutta 
Explicit 4thorder classical RungeKutta method 
expl_5_cash_karp_adptv 
Explicit adaptive 5thorder CashKarp method 
impl_1_euler_direct 
Implicit 1storder direct Euler method 
impl_1_euler_gmres 
Implicit 1storder Euler method with GMRES solver 
impl_2_midpoint_direct 
Implicit 2ndorder direct midpoint method 