#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Numerical solver of ordinary differential equations.
Solves the initial value problem for systems of first order
ordinary differential equations.
:Date: 2015-09-21
.. module:: ode
:platform: *nix, Windows
:synopsis: Numerical solver.
.. moduleauthor:: Daniel Weschke <daniel.weschke@directbox.de>
"""
from __future__ import division, print_function
from numpy import array, isnan, sum, zeros, dot
from numpy.linalg import norm, inv
[docs]def e1(f, x0, t, *p, verbose=False):
r"""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
:param f: the function to solve
:type f: function
:param x0: initial condition
:type x0: list
:param t: time
:type t: list
:param `*p`: parameters of the function (thickness, diameter,
...)
:param verbose: print information (default = False)
:type verbose: bool
Approximate the solution of the initial value problem
.. math ::
\dot{x} &= f(t,x) \\
x(t_0) &= x_0
Choose a value h for the size of every step and set
.. math ::
t_i = t_0 + i h ~,\quad i=1,2,\ldots,n
The derivative of the solution is approximated as the forward
difference equation
.. math ::
\dot{x}_i = f(t_i, x_i) = \frac{x_{i+1} - x_i}{t_{i+1}-t_i}
Therefore one step :math:`h` of the Euler method from
:math:`t_i` to :math:`t_{i+1}` is
.. math ::
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) \\
Example 1:
.. math ::
m\ddot{u} + d\dot{u} + ku = f(t) \\
\ddot{u} = m^{-1}(f(t) - d\dot{u} - ku) \\
with
.. math ::
x_1 &= u &\quad \dot{x}_1 = \dot{u} = x_2 \\
x_2 &= \dot{u} &\quad \dot{x}_2 = \ddot{u} \\
becomes
.. math ::
\dot{x}_1 &= x_2 \\
\dot{x}_2 &= m^{-1}(f(t) - d x_2 - k x_1) \\
or
.. math ::
\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}
Example 2:
.. math ::
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) \\
with
.. math ::
x_1 &= u &\quad \dot{x}_1 = \dot{u} = x_2 \\
x_2 &= \dot{u} &\quad \dot{x}_2 = \ddot{u} \\
becomes
.. math ::
\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) \\
or
.. math ::
\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}
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.
"""
x = zeros((len(t), len(x0))) # Preallocate array
x[0,:] = x0 # Initial condition gives solution at first t
for i in range(len(t)-1): # Calculation loop
Dt = t[i+1]-t[i]
dxdt = array(f(x[i,:], t[i], *p))
# Approximate solution at next value of x
x[i+1,:] = x[i,:] + dxdt*Dt
if verbose:
print('Numerical integration of ODE using explicit ' +
'first-order method (Euler / Runge-Kutta) was successful.')
return x
[docs]def e2(f, x0, t, *p, verbose=False):
r"""Explicit second-order method / Runge-Kutta 2nd order method.
:param f: the function to solve
:type f: function
:param x0: initial condition
:type x0: list
:param t: time
:type t: list
:param `*p`: parameters of the function (thickness, diameter,
...)
:param verbose: print information (default = False)
:type verbose: bool
"""
x = zeros((len(t), len(x0))) # Preallocate array
x[0,:] = x0 # Initial condition gives solution at first t
for i in range(len(t)-1): # Calculation loop
Dt = t[i+1]-t[i]
k_1 = array(f(x[i,:], t[i], *p))
k_2 = array(f(x[i,:]+0.5*Dt*k_1, t[i]+0.5*Dt, *p))
# Approximate solution at next value of x
x[i+1,:] = x[i,:] + k_2*Dt
if verbose:
print('Numerical integration of ODE using explicit ' +
'2th-order method (Runge-Kutta) was successful.')
return x
[docs]def e4(f, x0, t, *p, verbose=False):
r"""Explicit fourth-order method / Runge-Kutta 4th order method.
:param f: the function to solve
:type f: function
:param x0: initial condition
:type x0: list
:param t: time
:type t: list
:param `*p`: parameters of the function (thickness, diameter,
...)
:param verbose: print information (default = False)
:type verbose: bool
"""
x = zeros((len(t), len(x0))) # Preallocate array
x[0,:] = x0 # Initial condition
for i in range(len(t)-1): # Calculation loop
Dt = t[i+1]-t[i]
k_1 = array(f(x[i,:], t[i], *p))
k_2 = array(f(x[i,:]+0.5*Dt*k_1, t[i]+0.5*Dt, *p))
k_3 = array(f(x[i,:]+0.5*Dt*k_2, t[i]+0.5*Dt, *p))
k_4 = array(f(x[i,:]+k_3*Dt, t[i]+Dt, *p))
# Approximate solution at next value of x
x[i+1,:] = x[i,:] + 1./6*(k_1+2*k_2+2*k_3+k_4)*Dt
if verbose:
print('Numerical integration of ODE using explicit ' +
'4th-order method (Runge-Kutta) was successful.')
return x
[docs]def fpi(f, xi, ti, ti1, *p, max_iterations=1000, tol=1e-9,
verbose=False):
r"""Fixed-point iteration.
:param f: the function to iterate :math:`f = \dot{x}(x,t)`
:type f: function
:param xi: initial condition :math:`x_i`
:type xi: list
:param ti: time :math:`t_i`
:type ti: float
:param ti1: time :math:`t_{i+1}`
:type ti1: float
:param `*p`: parameters of the function (thickness, diameter,
...)
:param max_iterations: maximum number of iterations
:type max_iterations: int
:param tol: tolerance against residuum :math:`\varepsilon`
(default = 1e-9)
:type tol: float
:param verbose: print information (default = False)
:type verbose: bool
:returns: :math:`x_{i}`
.. math ::
x_{i,j=0} = x_{i}
.. math ::
x_{i,j+1} = x_i + \dot{x}(x_{i,j}, t_{i+1})\cdot(t_{i+1}-t_i)
.. math ::
\text{residuum} = \frac{\lVert x_{i,j+1}-x_{i,j}\rVert}
{\lVert x_{i,j+1} \rVert} < \varepsilon
.. math ::
x_{i} = x_{i,j=\text{end}}
"""
xij = xi
for j in range(max_iterations):
dxdt = array(f(xij, ti1, *p))
# Approximate solution at next value of x
xij1 = xi + dxdt * (ti1-ti)
residuum = norm(xij1-xij)/norm(xij1)
xij = xij1
if residuum < tol:
break
iterations = j+1 # number beginning with 1 therefore + 1
return xij, iterations
[docs]def i1(f, x0, t, *p, max_iterations=1000, tol=1e-9,
verbose=False):
r"""Implicite first-order method / backward Euler method.
:param f: the function to solve
:type f: function
:param x0: initial condition
:type x0: list
:param t: time
:type t: list
:param `*p`: parameters of the function (thickness, diameter,
...)
:param max_iterations: maximum number of iterations
:type max_iterations: int
:param tol: tolerance against residuum (default = 1e-9)
:type tol: float
:param verbose: print information (default = False)
:type verbose: bool
The backward Euler method has order one and is A-stable.
"""
iterations = zeros((len(t), 1))
x = zeros((len(t), len(x0))) # Preallocate array
x[0,:] = x0 # Initial condition gives solution at first t
# x(i+1) = x(i) + f(x(i+1), t(i+1)), exact value of
# f(x(i+1), t(i+1)) is not available therefore using
# Newton-Raphson method
for i in range(len(t)-1):
Dt = t[i+1]-t[i]
xi = x[i,:]
xi, iteration = fpi(f, xi, t[i], t[i+1], *p, max_iterations,
tol, verbose)
x[i+1,:] = xi
iterations[i] = iteration
if verbose:
print('Numerical integration of ODE using implicite ' +
'first-order method (Euler) was successful.')
return x, iterations
[docs]def newmark_newtonraphson(f, x0, xp0, xpp0, t, *p, gamma=.5,
beta=.25, max_iterations=1000, tol=1e-9, verbose=False):
r"""Newmark method.
:param f: the function to solve
:type f: function
:param x0: initial condition
:type x0: list
:param xp0: initial condition
:type xp0: list
:param xpp0: initial condition
:type xpp0: list
:param t: time
:type t: list
:param `*p`: parameters of the function (thickness, diameter,
...)
:param gamma: newmark parameter for velocity (default = 0.5)
:type gamma: float
:param beta: newmark parameter for displacement (default = 0.25)
:type beta: float
:param max_iterations: maximum number of iterations
:type max_iterations: int
:param tol: tolerance against residuum (default = 1e-9)
:type tol: float
:param verbose: print information (default = False)
:type verbose: bool
"""
iterations = zeros((len(t), 1))
x = zeros((len(t), len(x0))) # Preallocate array
xp = zeros((len(t), len(xp0))) # Preallocate array
xpp = zeros((len(t), len(xpp0))) # Preallocate array
x[0,:] = x0 # Initial condition gives solution at first t
xp[0,:] = xp0 # Initial condition gives solution at first t
xpp[0,:] = xpp0 # Initial condition gives solution at first t
for i in range(len(t)-1):
Dt = t[i+1]-t[i]
xi = x[i,:].reshape(3,1)
xpi = xp[i,:].reshape(3,1)
xppi = xpp[i,:].reshape(3,1)
x1 = xi
xp1 = xpi
xpp1 = xppi
j = 0
for j in range(max_iterations): # Fixed-point iteration
#dxdt = array(f(t[i+1], x1, p))
# Approximate solution at next value of x
#x11 = x[i,:] + dxdt*Dt
N, dN, dNp, dNpp = f(x1.reshape(-1,).tolist(),
xp1.reshape(-1,).tolist(), xpp1.reshape(-1,).tolist(),
t[i], *p)
if isnan(sum(dN)) or isnan(sum(dNp)) or isnan(sum(dNpp)):
print('divergiert')
break
xpp11 = xpp1 - dot(inv(dNpp), (N + dot(dN, (x1-xi)) + \
dot(dNp, (xp1-xpi))))
xp1 = xpi + Dt*( (1-gamma)*xppi + gamma*xpp11 )
x1 = xi + Dt*xpi + Dt**2*( (.5-beta)*xppi + beta*xpp11 )
residuum = norm(xpp11-xpp1)/norm(xpp11)
xpp1 = xpp11
if residuum < tol:
break
iterations[i] = j+1
xpp[i+1,:] = xpp1.reshape(-1,).tolist()
xp[i+1,:] = xp1.reshape(-1,).tolist()
x[i+1,:] = x1.reshape(-1,).tolist()
if verbose:
print('Numerical integration of ODE using explicite ' +
'newmark method was successful.')
return x, xp, xpp, iterations
# x = concatenate((x, xp, xpp), axis=1)
[docs]def newmark_newtonraphson_rdk(fnm, x0, xp0, xpp0, t, *p, gamma=.5,
beta=.25, max_iterations=1000, tol=1e-9, verbose=False):
r"""Newmark method.
:param f: the function to solve
:type f: function
:param x0: initial condition
:type x0: list
:param xp0: initial condition
:type xp0: list
:param xpp0: initial condition
:type xpp0: list
:param t: time
:type t: list
:param `*p`: parameters of the function (thickness, diameter,
...)
:param gamma: newmark parameter for velocity (default = 0.5)
:type gamma: float
:param beta: newmark parameter for displacement (default = 0.25)
:type beta: float
:param max_iterations: maximum number of iterations
:type max_iterations: int
:param tol: tolerance against residuum (default = 1e-9)
:type tol: float
:param verbose: print information (default = False)
:type verbose: bool
"""
iterations = zeros((len(t), 1))
x = zeros((len(t), len(x0))) # Preallocate array
xp = zeros((len(t), len(xp0))) # Preallocate array
xpp = zeros((len(t), len(xpp0))) # Preallocate array
x[0,:] = x0 # Initial condition gives solution at first t
xp[0,:] = xp0 # Initial condition gives solution at first t
xpp[0,:] = xpp0 # Initial condition gives solution at first t
for i in range(len(t)-1):
Dt = t[i+1]-t[i]
rm, rmx, rmxpp, rd, rdx, rdxp, rk, rkx, f = fnm(x[i,:],
xp[i,:], xpp[i,:], t[i], *p)
xi = x[i,:].reshape(3,1)
xpi = xp[i,:].reshape(3,1)
xppi = xpp[i,:].reshape(3,1)
x1 = xi
xp1 = xpi
xpp1 = xppi
j = 0
for j in range(max_iterations): # Fixed-point iteration
#dxdt = array(f(t[i+1], x1, p))
# Approximate solution at next value of x
#x11 = x[i,:] + dxdt*Dt
r = (rmx+rdx+rkx)*Dt**2./4 + rdxp*Dt/2 + rmxpp
rp = f - (rm + dot(rmx, (Dt*xpi+Dt**2./4*xppi)) - \
dot(rmxpp, xppi) + \
rd + dot(rdx, (Dt*xpi+Dt**2./4*xppi)) + \
dot(rdxp, Dt/2*xppi) + \
rk + dot(rkx, (Dt*xpi+Dt**2./4*xppi)) )
xpp11 = dot(inv(r), rp)
xp1 = xpi + Dt*( (1-gamma)*xppi + gamma*xpp11 )
x1 = xi + Dt*xpi + Dt**2*( (.5-beta)*xppi + beta*xpp11 )
residuum = norm(xpp11-xpp1)/norm(xpp11)
xpp1 = xpp11
if residuum < tol:
break
iterations[i] = j+1
xpp[i+1,:] = xpp1.reshape(-1,).tolist()
xp[i+1,:] = xp1.reshape(-1,).tolist()
x[i+1,:] = x1.reshape(-1,).tolist()
if verbose:
print('Numerical integration of ODE using explicite ' +
'newmark method was successful.')
return x, xp, xpp, iterations
# x = concatenate((x, xp, xpp), axis=1)