+#!/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)
+
