pylib.numerical.ode module

Numerical solver of ordinary differential equations.

Solves the initial value problem for systems of first order ordinary differential equations.

Date

2015-09-21

Approximate the solution \(x(t)\) of the initial value problem

\[\begin{split}\frac{\mathrm{d}x}{\mathrm{d}t} = \dot{x} &= f(t,x) \\ x(t_0) &= x_0 \\ t &\in [t_0, t_n]\end{split}\]

e1(f, x0, t, *p, verbose=False)

Explicit first-order method / (standard, or forward) Euler method / Runge-Kutta 1st order method.

de: Euler’sche Polygonzugverfahren / explizite Euler-Verfahren / Euler-Cauchy-Verfahren / Euler-vorwärts-Verfahren

Parameters

• f (function) – the function to solve

• x0 (list) – initial condition

• t (list) – time

• *p – parameters of the function (thickness, diameter, …)

• verbose (bool) – print information (default = False)

Approximate the solution \(x(t)\) of the initial value problem

\[\begin{split}\dot{x} &= f(t,x) \\ x(t_0) &= x_0 \\ t &\in [t_0, t_n]\end{split}\]

Choose a value h for the size of every step and set

\[t_{i+1} = t_0 + i h = t_i + h ~,\quad i=0,1,2,\ldots,n-1\]

The derivative of the solution is approximated as the forward difference equation

\[\dot{x}_i = f(t_i, x_i) = \frac{x_{i+1} - x_i}{t_{i+1}-t_i}\]

Therefore one step \(h\) of the Euler method from \(t_i\) to \(t_{i+1}\) is

\[\begin{split}x_{i+1} &= x_i + (t_{i+1}-t_i) f(t_i, x_i) \\ x_{i+1} &= x_i + h f(t_i, x_i) \\\end{split}\]

Example 1:

\[\begin{split}m\ddot{u} + d\dot{u} + ku = f(t) \\ \ddot{u} = m^{-1}(f(t) - d\dot{u} - ku) \\\end{split}\]

with

\[\begin{split}x_1 &= u &\quad \dot{x}_1 = \dot{u} = x_2 \\ x_2 &= \dot{u} &\quad \dot{x}_2 = \ddot{u} \\\end{split}\]

becomes

\[\begin{split}\dot{x}_1 &= x_2 \\ \dot{x}_2 &= m^{-1}(f(t) - d x_2 - k x_1) \\\end{split}\]

or

\[\begin{split}\dot{x} &= f(t,x) \\ \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} &= \begin{bmatrix} x_2 \\ m^{-1}(f(t) - d x_2 - k x_1) \end{bmatrix} \\ &= \begin{bmatrix} 0 \\ m^{-1} f(t) \end{bmatrix} + \begin{bmatrix} 0 & 1 \\ -m^{-1} k & -m^{-1} d \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}\end{split}\]

Example 2:

\[\begin{split}m(u)\ddot{u} + d(u,\dot{u})\dot{u} + k(u)u = f(t) \\ \ddot{u} = m^{-1}(u)(f(t) - d(u,\dot{u})\dot{u} - k(u)u) \\\end{split}\]

with

\[\begin{split}x_1 &= u &\quad \dot{x}_1 = \dot{u} = x_2 \\ x_2 &= \dot{u} &\quad \dot{x}_2 = \ddot{u} \\\end{split}\]

becomes

\[\begin{split}\dot{x}_1 &= x_2 \\ \dot{x}_2 &= m^{-1}(x_1)(f(t) - d(x_1,x_2) x_2 - k(x_1) x_1) \\\end{split}\]

or

\[\begin{split}\dot{x} &= f(t,x) \\ \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} &= \begin{bmatrix} x_2 \\ m^{-1}(x_1)(f(t) - d(x_1,x_2) x_2 - k(x_1) x_1) \end{bmatrix} \\ &= \begin{bmatrix} 0 \\ m^{-1}(x_1) f(t) \end{bmatrix} + \begin{bmatrix} 0 & 1 \\ -m^{-1}(x_1) k(x_1) & -m^{-1} d(x_1,x_2) \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}\end{split}\]

The Euler method is a first-order method, which means that the local error (error per step) is proportional to the square of the step size, and the global error (error at a given time) is proportional to the step size.

e2(f, x0, t, *p, verbose=False)

Explicit second-order method / Runge-Kutta 2nd order method.

Parameters

• f (function) – the function to solve

• x0 (list) – initial condition

• t (list) – time

• *p – parameters of the function (thickness, diameter, …)

• verbose (bool) – print information (default = False)

e4(f, x0, t, *p, verbose=False)

Explicit fourth-order method / Runge-Kutta 4th order method.

Parameters

• f (function) – the function to solve

• x0 (list) – initial condition

• t (list) – time

• *p – parameters of the function (thickness, diameter, …)

• verbose (bool) – print information (default = False)

Problem

\[\begin{split}\dot{y} &= f(t, y) \\ y(t_0) &= y_0 \\ t &\in [t_0, t_n]\end{split}\]

Increment \(\delta t = t_{i+1}-t_i ~,~~ i=0,1,2,\ldots,n-1\)

\[\begin{split}y_{n+1} &= y_{i} + \tfrac{1}{6}( \delta y_{i,1} + 2\delta y_{i,2} + 2\delta y_{i,3} + \delta y_{i,4}) \\ & \qquad \text{with} \\ \delta y_{i,1} &= \delta t \cdot y'(t_{i}, ~ y_{i}) \\ \delta y_{i,2} &= \delta t \cdot y'(t_{i}+\tfrac{1}{2}\delta t, ~ y_{i}+\tfrac{1}{2}\delta y_{i,1}) \\ \delta y_{i,3} &= \delta t \cdot y'(t_{i}+\tfrac{1}{2}\delta t, ~ y_{i}+\tfrac{1}{2}\delta y_{i,2}) \\ \delta y_{i,4} &= \delta t \cdot y'(t_{i}+\delta t, ~ y_{i}+\delta y_{i,3})\end{split}\]

fpi(f, xi, ti, ti1, *p, max_iterations=1000, tol=1e-09, verbose=False)

Fixed-point iteration.

Parameters

• f (function) – the function to iterate \(f = \dot{x}(x,t)\)

• xi (list) – initial condition \(x_i\)

• ti (float) – time \(t_i\)

• ti1 (float) – time \(t_{i+1}\)

• *p – parameters of the function (thickness, diameter, …)

• max_iterations (int) – maximum number of iterations

• tol (float) – tolerance against residuum \(\varepsilon\) (default = 1e-9)

• verbose (bool) – print information (default = False)

Returns

\(x_{i}\)

\[x_{i,j=0} = x_{i}\]

\[x_{i,j+1} = x_i + \dot{x}(x_{i,j}, t_{i+1})\cdot(t_{i+1}-t_i)\]

\[\text{residuum} = \frac{\lVert x_{i,j+1}-x_{i,j}\rVert} {\lVert x_{i,j+1} \rVert} < \varepsilon\]

\[x_{i} = x_{i,j=\text{end}}\]

i1(f, x0, t, *p, max_iterations=1000, tol=1e-09, verbose=False)

Implicite first-order method / backward Euler method.

Parameters

• f (function) – the function to solve

• x0 (list) – initial condition

• t (list) – time

• *p – parameters of the function (thickness, diameter, …)

• max_iterations (int) – maximum number of iterations

• tol (float) – tolerance against residuum (default = 1e-9)

• verbose (bool) – print information (default = False)

The backward Euler method has order one and is A-stable.

newmark_newtonraphson(f, x0, xp0, xpp0, t, *p, gamma=0.5, beta=0.25, max_iterations=1000, tol=1e-09, verbose=False)

Newmark method.

Parameters

• f (function) – the function to solve

• x0 (list) – initial condition

• xp0 (list) – initial condition

• xpp0 (list) – initial condition

• t (list) – time

• *p – parameters of the function (thickness, diameter, …)

• gamma (float) – newmark parameter for velocity (default = 0.5)

• beta (float) – newmark parameter for displacement (default = 0.25)

• max_iterations (int) – maximum number of iterations

• tol (float) – tolerance against residuum (default = 1e-9)

• verbose (bool) – print information (default = False)

newmark_newtonraphson_mdk(fmdk, x0, xp0, xpp0, t, *p, gamma=0.5, beta=0.25, max_iterations=1000, tol=1e-09, verbose=False)

Newmark method.

Using m mass, d damping and k stiffness formulation.

Parameters

• f (function) – the function to solve

• x0 (list) – initial condition

• xp0 (list) – initial condition

• xpp0 (list) – initial condition

• t (list) – time

• *p – parameters of the function (thickness, diameter, …)

• gamma (float) – newmark parameter for velocity (default = 0.5)

• beta (float) – newmark parameter for displacement (default = 0.25)

• max_iterations (int) – maximum number of iterations

• tol (float) – tolerance against residuum (default = 1e-9)

• verbose (bool) – print information (default = False)