@@ -175,6 +210,12 @@
S
geometry module
-
2D
-
@date: 2019-03-21
-@author: Daniel Weschke
-
+2D geometry objects.
+
+Date
+2019-03-21
+
+
cubic( point1 , angle1 , point2 , angle2 , samples=50 )
-returns
-([sample_point1_x, sample_point2_x, …], [sample_points1_y, sample_point2_y, …])
+
+Parameters
+
+point1 (tuple ) – one end point
angle1 (int or float ) – the slope at the one end point
point2 (tuple ) – other end point
angle2 (int or float ) – the slope at the other end point
samples (int ) – number of sampling points
+Returns
+([sample_point1_x, sample_point2_x, …],
+[sample_points1_y, sample_point2_y, …])
+Return type
+tuple
+
cubic_deg( point1 , angle1 , point2 , angle2 )
-returns
-([sample_point1_x, sample_point2_x, …], [sample_points1_y, sample_point2_y, …])
+
+Parameters
+
+point1 (tuple ) – one end point
angle1 (int or float ) – the slope at the one end point
point2 (tuple ) – other end point
angle2 (int or float ) – the slope at the other end point
+Returns
+([sample_point1_x, sample_point2_x, …],
+[sample_points1_y, sample_point2_y, …])
+Return type
+tuple
+
+
line( point1 , point2 , samples=2 )
-samples: number of sampling points
-y = (y2-y1)/(x2-x1)*(x-x1) + y1
-
-returns ((point1_x, point2_x), (points1_y, point2_y))
-or
-([sample_point1_x, sample_point2_x, …], [sample_points1_y, sample_point2_y, …])
+
+\[y = \frac{y_2-y_1}{x_2-x_1}(x-x_1) + y_1\]
+
+Parameters
+
+point1 (tuple ) – one end point
point2 (tuple ) – other end point
samples (int ) – number of sampling points
+Returns
+((point1_x, point2_x), (points1_y, point2_y)) or
+([sample_point1_x, sample_point2_x, …],
+[sample_points1_y, sample_point2_y, …])
+Return type
+tuple
@@ -66,26 +116,57 @@ or
points( *pts )
-returns
-((point1_x, point2_x), (point1_y, point2_y), …)
+
+Parameters
+*pts – points to rearrange
+Returns
+((point1_x, point2_x), (point1_y, point2_y), …)
+Return type
+tuple
+
rectangle( width , height )
-returns
-(point1, point2, point3, point4)
+
+Parameters
+
+
+Returns
+(point1, point2, point3, point4)
+Return type
+tuple
+
rotate( origin , angle , *pts , **kwargs )
-Rotate a point or polygon counterclockwise by a given angle around a given origin.
-The angle should be given in radians.
-
-returns (point_x, point_y)
-or
-(point1, point2, …)
+Rotate a point or polygon counterclockwise by a given angle around a given
+origin. The angle should be given in radians.
+
+Parameters
+
+
+Returns
+(point_x, point_y) or (point1, point2, …)
+Return type
+tuple
@@ -93,31 +174,66 @@ or
rotate_deg( origin , angle , *pts , **kwargs )
-Rotate a point or polygon counterclockwise by a given angle around a given origin.
-The angle should be given in degrees.
-
-returns (point_x, point_y)
-or
-(point1, point2, …)
+Rotate a point or polygon counterclockwise by a given angle around a given
+origin. The angle should be given in degrees.
+
+Parameters
+
+
+Returns
+(point_x, point_y) or (point1, point2, …)
+Return type
+tuple
+
square( width )
-returns
-(point1, point2, point3, point4)
+
+Parameters
+width (int or float ) – the edge size of the square
+Returns
+(point1, point2, point3, point4)
+Return type
+tuple
+
+
translate( vec , *pts )
Translate a point or polygon by a given vector.
-
-returns (point_x, point_y)
-or
-(point1, point2, …)
+
+Parameters
+
+
+Returns
+(point_x, point_y) or (point1, point2, …)
+Return type
+tuple
diff --git a/docs/build/html/index.html b/docs/build/html/index.html
index cccc0f5..88d9565 100644
--- a/docs/build/html/index.html
+++ b/docs/build/html/index.html
@@ -13,6 +13,7 @@
+
date module fit module geometry module solver module solver_model module
diff --git a/docs/build/html/objects.inv b/docs/build/html/objects.inv
index 13d3524..ab7137e 100644
--- a/docs/build/html/objects.inv
+++ b/docs/build/html/objects.inv
@@ -2,9 +2,10 @@
# Project: pylib
# Version:
# The remainder of this file is compressed using zlib.
-xڕAN0E9lS7`Ăe4'ǑȎk0܄`I7Ul"{#kmmD3m w,)I$#E-.FFmmpQ)LY!4WTy@mQZ)[F]WcaB^ilzZ!)R~?kcytڸM S.6g$?Mx!Ztcv J(aGa
т]QYФ
-Cȉ0!?XqS
-G
-*< 4
-? <؊ш?Nn>
-D4swD71QݹɼoMQ)NCtKӴj 0YT5"ϟS- uX=0y4h
\ No newline at end of file
+xڥMr s
+RVS++ Tʋ,)Zb~%$AxPA^CZyFӝbP@m
+5H;Uȷ^ZCn2E:F[YͤF|73`0Pү.tCN
+vƵ5[p
ڨkӠz0 xu ϋ2XX wa +[0 X
ޝM.|8H>-UK
+%J{+ztjC
+g*lm
m 3jV}Wr
+-x 3sPʼ~h h8 њR' fL8)$gū*j? K;?A/ɾIJ_ {SMM'&#gNCjn&c ṯ㗧x3/x*J ~?>! ؠC#gBovfem&էt*Sj̳;tMܧϬ9wsL?
\ No newline at end of file
diff --git a/docs/build/html/py-modindex.html b/docs/build/html/py-modindex.html
index dcada82..180c2c8 100644
--- a/docs/build/html/py-modindex.html
+++ b/docs/build/html/py-modindex.html
@@ -13,6 +13,7 @@
+
@@ -54,29 +56,42 @@
- data
-
+ data(*nix, Windows)
+ Handle data files.
- date
-
+ date(*nix, Windows)
+ Special dates.
f
- fit
-
+ fit(*nix, Windows)
+ Function and approximation.
g
- geometry
-
+ geometry(*nix, Windows)
+ Geometry objects.
+
+ s
+
+ solver(*nix, Windows)
+ Numerical solver.
+
+
+ solver_model(*nix, Windows)
+ Mechanical models.
diff --git a/docs/build/html/search.html b/docs/build/html/search.html
index bd5a461..4784e07 100644
--- a/docs/build/html/search.html
+++ b/docs/build/html/search.html
@@ -14,6 +14,7 @@
+
solver module — pylib 2019.5.19 documentation
+
+
+
+
+
+
+
+
+
solver module
+
Numerical solver of ordinary differential equations.
+
Solves the initial value problem for systems of first order ordinary differential
+equations.
+
+Date
+2015-09-21
+
+
+
+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 of the initial value problem
+
+\[\begin{split}\dot{x} &= f(t,x) \\
+x(t_0) &= x_0\end{split}\]
+Choose a value h for the size of every step and set
+
+\[t_i = t_0 + i h ~,\quad i=1,2,\ldots,n\]
+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}_1 \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)
+
+
+
+
+
+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_rdk( fnm , x0 , xp0 , xpp0 , t , *p , gamma=0.5 , beta=0.25 , maxIterations=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)
+
+
+
+
+
+
+
+
+
+
+
+
+
solver_model module — pylib 2019.5.19 documentation
+
+
+
+
+
+
+
+
+
solver_model module
+
Mathmatical models governed by ordinary differential equations.
+
Describes initial value problems as systems of first order ordinary differential
+equations.
+
+Date
+2019-05-25
+
+
+
+disk( x , t , *p )
+Rotation of an eccentric disk.
+
+Parameters
+
+
+
+
+
+
+
+disk_nm( xn , xpn , xppn , t , *p )
+Rotation of an eccentric disk.
+
+Parameters
+
+xn (list ) – values of the function
xpn (list ) – first derivative values of the function
xppn (list ) – second derivative values of the function
t (list ) – time
*p –
parameters of the function
+
+diameter
eccentricity
torque
+
+
+
+
+
+
+disk_nmmdk( xn , xpn , xppn , t , *p )
+Rotation of an eccentric disk.
+
+Parameters
+
+xn (list ) – values of the function
xpn (list ) – derivative values of the function
xppn (list ) – second derivative values of the function
t (list ) – time
*p –
parameters of the function
+
+diameter
eccentricity
torque
+
+
+
+
+
+
+
+
+
+
+
+
+
+
src package — pylib 2019.5.19 documentation
-
-
-
-
-
-
-
-
-
src package
-
-
Submodules
-
-
-
src.data module
-
-
-
src.date module
-
-
-
src.fit module
-
-
-
src.geometry module
-
-
-
Module contents
-
-
-
-
-
-
-
-
-
-
-
+"""
from __future__ import print_function
import pickle
def data_read(file_name, x_column, y_column):
- """read ascii data file"""
+ """Read ascii data file.
+
+ :param filename: file to read
+ :type filename: str
+ :param x_column: column index for the x data (first column is 0)
+ :type x_column: int
+ :param y_column: column index for the y data (first column is 0)
+ :type y_column: int
+
+ :returns: x and y
+ :rtype: tuple(list, list)
+ """
import re
file = open(file_name)
x = []
@@ -15,7 +36,16 @@ def data_read(file_name, x_column, y_column):
return x, y
def data_load(file_name, verbose=False):
- """load stored program objects from binary file"""
+ """Load stored program objects from binary file.
+
+ :param file_name: file to load
+ :type file_name: str
+ :param verbose: verbose information (default = False)
+ :type verbose: bool
+
+ :returns: loaded data
+ :rtype: object
+ """
if verbose:
print('check if data is available')
try:
@@ -31,11 +61,18 @@ def data_load(file_name, verbose=False):
return object_data
def data_store(file_name, object_data):
- """store program objects to binary file"""
+ """Store program objects to binary file.
+
+ :param file_name: file to store
+ :type file_name: str
+ :param object_data: data to store
+ :type object_data: object
+ """
with open(file_name, 'wb') as output:
pickle.dump(object_data, output, pickle.HIGHEST_PROTOCOL) # every dump needs a load
def main():
+ """Main function."""
file_name = "slit_test_scan.dat"
x, y = data_read(file_name, 3, 2)
print(x)
diff --git a/src/date.py b/src/date.py
index ee7ad8a..476ee43 100644
--- a/src/date.py
+++ b/src/date.py
@@ -1,31 +1,51 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
-"""
-Created on Mon Jan 15 21:52:34 2018
+"""Calculate spacial dates.
-@author: Daniel Weschke
+:Date: 2018-01-15
+
+.. module:: date
+ :platform: *nix, Windows
+ :synopsis: Special dates.
+
+.. moduleauthor:: Daniel Weschke
"""
from __future__ import division, print_function, unicode_literals
def gaußsche_osterformel(year):
"""Gaußsche Osterformel.
+
+ :param year: the year to calculate the Easter Sunday
+ :type year: int
- :returns: int -- the day of Easter Sunday as a day in march.
+ :returns: the day of Easter Sunday as a day in march.
+ :rtype: int
+
+ :ivar X: Das Jahr / year
+ :vartype X: int
+ :ivar K(X): Die Säkularzahl
+ :vartype K(X): int
+ :ivar M(X): Die säkulare Mondschaltung
+ :vartype M(X): int
+ :ivar S(K): Die säkulare Sonnenschaltung
+ :vartype S(K): int
+ :ivar A(X): Den Mondparameter
+ :vartype A(X): int
+ :ivar D(A,M): Den Keim für den ersten Vollmond im Frühling
+ :vartype D(A,M): int
+ :ivar R(D,A): Die kalendarische Korrekturgröße
+ :vartype R(D,A): int
+ :ivar OG(D,R): Die Ostergrenze
+ :vartype OG(D,R): int
+ :ivar SZ(X,S): Den ersten Sonntag im März
+ :vartype SZ(X,S): int
+ :ivar OE(OG,SZ): Die Entfernung des Ostersonntags von der Ostergrenze (Osterentfernung in Tagen)
+ :vartype OE(OG,SZ): int
+ :ivar OS(OG,OE): Das Datum des Ostersonntags als Märzdatum (32. März = 1. April usw.)
+ :vartype OS(OG,OE): int
Algorithmus gilt für den Gregorianischen Kalender.
- * X -- Das Jahr / year
- * K(X) -- Die Säkularzahl
- * M(K) -- Die säkulare Mondschaltung
- * S(K) -- Die säkulare Sonnenschaltung
- * A(X) -- Den Mondparameter
- * D(A,M) -- Den Keim für den ersten Vollmond im Frühling
- * R(D,A) -- Die kalendarische Korrekturgröße
- * OG(D,R) -- Die Ostergrenze
- * SZ(X,S) -- Den ersten Sonntag im März
- * OE(OG,SZ) -- Die Entfernung des Ostersonntags von der Ostergrenze (Osterentfernung in Tagen)
- * OS(OG,OE) -- Das Datum des Ostersonntags als Märzdatum (32. März = 1. April usw.)
-
source: https://de.wikipedia.org/wiki/Gau%C3%9Fsche_Osterformel
"""
x = year
@@ -43,8 +63,12 @@ def gaußsche_osterformel(year):
def easter_sunday(year):
"""Easter Sunday.
+
+ :param year: the year to calculate the Easter Sunday
+ :type year: int
- :returns: datetime.date -- the day of Easter Sunday"""
+ :returns: the day of Easter Sunday
+ :rtype: datetime.date"""
import datetime
march = datetime.date(year, 3, 1)
day = march + datetime.timedelta(days=gaußsche_osterformel(year))
@@ -52,32 +76,48 @@ def easter_sunday(year):
def easter_friday(year):
"""Easter Friday.
+
+ :param year: the year to calculate the Easter Friday
+ :type year: int
- :returns: datetime.date -- the day of Easter Friday"""
+ :returns: the day of Easter Friday
+ :rtype: datetime.date"""
import datetime
day = easter_sunday(year) + datetime.timedelta(days=-2)
return day
def easter_monday(year):
"""Easter Monday.
+
+ :param year: the year to calculate the Easter Monday
+ :type year: int
- :returns: datetime.date -- the day of Easter Monday"""
+ :returns: the day of Easter Monday
+ :rtype: datetime.date"""
import datetime
day = easter_sunday(year) + datetime.timedelta(days=+1)
return day
def ascension_of_jesus(year):
"""Ascension of Jesus.
+
+ :param year: the year to calculate the ascension of Jesus
+ :type year: int
- :returns: datetime.date -- the day of ascension of Jesus"""
+ :returns: the day of ascension of Jesus
+ :rtype: datetime.date"""
import datetime
day = easter_sunday(year) + datetime.timedelta(days=+39)
return day
def pentecost(year):
"""Pentecost.
+
+ :param year: the year to calculate the Pentecost
+ :type year: int
- :returns: datetime.date -- the day of Pentecost"""
+ :returns: the day of Pentecost
+ :rtype: datetime.date"""
import datetime
day = easter_sunday(year) + datetime.timedelta(days=+49)
return day
diff --git a/src/fit.py b/src/fit.py
index fef8b48..781dc5d 100644
--- a/src/fit.py
+++ b/src/fit.py
@@ -1,3 +1,13 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Function and approximation.
+
+.. module:: fit
+ :platform: *nix, Windows
+ :synopsis: Function and approximation.
+
+.. moduleauthor:: Daniel Weschke
+"""
from __future__ import print_function
from pylab import array, argmax, subplot, plot, title, xlim, show, gradient, exp, sqrt, log, linspace
from scipy.optimize import curve_fit
@@ -6,15 +16,21 @@ from data import *
def gauss(x, *p):
"""Gauss distribution function.
+ .. math::
+ f(x)=ae^{-(x-b)^{2}/(2c^{2})}
+
:param x: positions where the gauss function will be calculated
+ :type x: int or float or list or numpy.ndarray
:param p: gauss parameters [a, b, c, d]:
- * a -- amplitude (integral = 1 if a = 1/(c*sqrt(2*pi)))
- * b -- expected value mu (position of maximum, default = 0)
- * c -- standard deviation sigma (variance sigma**2 = c**2)
+ * a -- amplitude (:math:`\int y \\,\\mathrm{d}x=1 \Leftrightarrow a=1/(c\\sqrt{2\\pi})` )
+ * b -- expected value :math:`\\mu` (position of maximum, default = 0)
+ * c -- standard deviation :math:`\\sigma` (variance :math:`\\sigma^2=c^2`)
* d -- vertical offset (default = 0)
+ :type p: list
- :returns: array -- gauss values at given positions x
+ :returns: gauss values at given positions x
+ :rtype: numpy.ndarray
"""
x = array(x) # cast e. g. list to numpy array
a, b, c, d = p
@@ -24,15 +40,27 @@ def gauss_fit(x, y, e=None, x_fit=None, verbose=False):
"""Fit Gauss distribution function to data.
:param x: positions
+ :type x: int or float or list or numpy.ndarray
:param y: values
- :param e: error values
- :param x_fit: positions of fitted function (default steps is 3*len(x) but min 150)
+ :type y: int or float or list or numpy.ndarray
+ :param e: error values (default = None)
+ :type e: int or float or list or numpy.ndarray
+ :param x_fit: positions of fitted function (default = None, if None then x
+ is used)
+ :type x_fit: int or float or list or numpy.ndarray
+ :param verbose: verbose information (default = False)
+ :type verbose: bool
:returns:
- * y_fit -- values
- * popt -- parameters of gauss distribution function (amplitude a, expected
- value mu, standard deviation sigma, vertical offset d)
- * FWHM -- full width at half maximum
+ * numpy.ndarray -- fitted values (y_fit)
+ * numpy.ndarray -- parameters of gauss distribution function (popt:
+ amplitude a, expected value :math:`\\mu`, standard deviation
+ :math:`\\sigma`, vertical offset d)
+ * numpy.float64 -- full width at half maximum (FWHM)
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`gauss`
"""
x = array(x) # cast e. g. list to numpy array
y = array(y) # cast e. g. list to numpy array
diff --git a/src/geometry.py b/src/geometry.py
index b0f19a1..d3099ce 100644
--- a/src/geometry.py
+++ b/src/geometry.py
@@ -1,37 +1,46 @@
+#!/usr/bin/env python
# -*- coding: utf-8 -*-
-"""
-2D
+"""2D geometry objects.
-@date: 2019-03-21
-@author: Daniel Weschke
+:Date: 2019-03-21
+
+.. module:: geometry
+ :platform: *nix, Windows
+ :synopsis: Geometry objects.
+
+.. moduleauthor:: Daniel Weschke
"""
import math
import numpy as np
def translate(vec, *pts):
- """\
- Translate a point or polygon by a given vector.
+ """Translate a point or polygon by a given vector.
- returns
- (point_x, point_y)
- or
- (point1, point2, ...)
+ :param vec: translation vector
+ :type vec: tuple
+ :param `*pts`: points to translate
+
+ :returns: (point_x, point_y) or (point1, point2, ...)
+ :rtype: tuple
"""
vx, vy = vec
return tuple([(x+vx, y+vy) for (x, y) in pts])
def rotate(origin, angle, *pts, **kwargs):
- """\
- Rotate a point or polygon counterclockwise by a given angle around a given origin.
+ """Rotate a point or polygon counterclockwise by a given angle around a given
+ origin. The angle should be given in radians.
- The angle should be given in radians.
+ :param origin: the center of rotation
+ :type origin: tuple
+ :param angle: the rotation angle
+ :type angle: int or float
+ :param `*pts`: points to rotate
+ :param `**kwargs`: options
- returns
- (point_x, point_y)
- or
- (point1, point2, ...)
+ :returns: (point_x, point_y) or (point1, point2, ...)
+ :rtype: tuple
"""
ox, oy = origin
@@ -50,23 +59,34 @@ def rotate(origin, angle, *pts, **kwargs):
def rotate_deg(origin, angle, *pts, **kwargs):
- """\
- Rotate a point or polygon counterclockwise by a given angle around a given origin.
+ """Rotate a point or polygon counterclockwise by a given angle around a given
+ origin. The angle should be given in degrees.
- The angle should be given in degrees.
+ :param origin: the center of rotation
+ :type origin: tuple
+ :param angle: the rotation angle
+ :type angle: int or float
+ :param `*pts`: points to rotate
+ :param `**kwargs`: options
- returns
- (point_x, point_y)
- or
- (point1, point2, ...)
+ :returns: (point_x, point_y) or (point1, point2, ...)
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`rotate`
"""
return rotate(origin, angle*math.pi/180, *pts, **kwargs)
def rectangle(width, height):
"""\
- returns
- (point1, point2, point3, point4)
+ :param width: the width of the rectangle
+ :type width: int or float
+ :param height: the height of the rectangle
+ :type height: int or float
+
+ :returns: (point1, point2, point3, point4)
+ :rtype: tuple
"""
pt1 = (-width/2, -height/2)
pt2 = (width/2, -height/2)
@@ -77,8 +97,14 @@ def rectangle(width, height):
def square(width):
"""\
- returns
- (point1, point2, point3, point4)
+ :param width: the edge size of the square
+ :type width: int or float
+
+ :returns: (point1, point2, point3, point4)
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`rectangle`
"""
return rectangle(width, width)
@@ -89,22 +115,30 @@ def square(width):
def points(*pts):
"""\
- returns
- ((point1_x, point2_x), (point1_y, point2_y), ...)
+ :param `*pts`: points to rearrange
+
+ :returns: ((point1_x, point2_x), (point1_y, point2_y), ...)
+ :rtype: tuple
"""
return zip(*pts)
def line(point1, point2, samples=2):
"""\
- samples: number of sampling points
+ .. math::
+ y = \\frac{y_2-y_1}{x_2-x_1}(x-x_1) + y_1
+
+ :param point1: one end point
+ :type point1: tuple
+ :param point2: other end point
+ :type point2: tuple
+ :param samples: number of sampling points
+ :type samples: int
- y = (y2-y1)/(x2-x1)*(x-x1) + y1
-
- returns
- ((point1_x, point2_x), (points1_y, point2_y))
- or
- ([sample_point1_x, sample_point2_x, ...], [sample_points1_y, sample_point2_y, ...])
+ :returns: ((point1_x, point2_x), (points1_y, point2_y)) or
+ ([sample_point1_x, sample_point2_x, ...],
+ [sample_points1_y, sample_point2_y, ...])
+ :rtype: tuple
"""
p1x, p1y = point1
p2x, p2y = point2
@@ -122,8 +156,20 @@ def line(point1, point2, samples=2):
def cubic(point1, angle1, point2, angle2, samples=50):
"""\
- returns
- ([sample_point1_x, sample_point2_x, ...], [sample_points1_y, sample_point2_y, ...])
+ :param point1: one end point
+ :type point1: tuple
+ :param angle1: the slope at the one end point
+ :type angle1: int or float
+ :param point2: other end point
+ :type point2: tuple
+ :param angle2: the slope at the other end point
+ :type angle2: int or float
+ :param samples: number of sampling points
+ :type samples: int
+
+ :returns: ([sample_point1_x, sample_point2_x, ...],
+ [sample_points1_y, sample_point2_y, ...])
+ :rtype: tuple
"""
p1x, p1y = point1
p2x, p2y = point2
@@ -142,7 +188,20 @@ def cubic(point1, angle1, point2, angle2, samples=50):
def cubic_deg(point1, angle1, point2, angle2):
"""\
- returns
- ([sample_point1_x, sample_point2_x, ...], [sample_points1_y, sample_point2_y, ...])
+ :param point1: one end point
+ :type point1: tuple
+ :param angle1: the slope at the one end point
+ :type angle1: int or float
+ :param point2: other end point
+ :type point2: tuple
+ :param angle2: the slope at the other end point
+ :type angle2: int or float
+
+ :returns: ([sample_point1_x, sample_point2_x, ...],
+ [sample_points1_y, sample_point2_y, ...])
+ :rtype: tuple
+
+ .. seealso::
+ :meth:`cubic`
"""
return cubic(point1, angle1 * math.pi/180, point2, angle2 * math.pi/180)
diff --git a/src/solver.py b/src/solver.py
new file mode 100644
index 0000000..0e5fc10
--- /dev/null
+++ b/src/solver.py
@@ -0,0 +1,357 @@
+#!/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:: solver
+ :platform: *nix, Windows
+ :synopsis: Numerical solver.
+
+.. moduleauthor:: Daniel Weschke
+"""
+
+from __future__ import division, print_function
+from numpy import array, isnan, sum, zeros, dot
+from numpy.linalg import norm, inv
+
+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
+
+ 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}_1 \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 t=0.
+ for i in range(len(t)-1):
+ Dt = t[i+1]-t[i]
+ dxdt = array(f(x[i,:], t[i], *p))
+ x[i+1,:] = x[i,:] + dxdt*Dt # Approximate solution at next value of x
+ if verbose:
+ print('Numerical integration of ODE using explicit first-order method (Euler / Runge-Kutta) was successful.')
+ return x
+
+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)))
+ 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))
+ x[i+1,:] = x[i,:] + k_2*Dt # main equation
+ if verbose:
+ print('Numerical integration of ODE using explicit 2th-order method (Runge-Kutta) was successful.')
+ return x
+
+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 gives solution at t=0.
+ 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))
+ x[i+1,:] = x[i,:] + 1./6*(k_1+2*k_2+2*k_3+k_4)*Dt # main equation
+ if verbose:
+ print('Numerical integration of ODE using explicit 4th-order method (Runge-Kutta) was successful.')
+ return x
+
+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 t=0.
+ for i in range(len(t)-1):
+ Dt = t[i+1]-t[i]
+ x1 = x[i,:]
+ for j in range(max_iterations): # fixed-point iteration
+ dxdt = array(f(x1, t[i+1], *p))
+ x11 = x[i,:] + dxdt*Dt # Approximate solution at next value of x
+ residuum = norm(x11-x1)/norm(x11)
+ x1 = x11
+ if residuum < tol:
+ break
+ iterations[i] = j+1
+ x[i+1,:] = x1
+ if verbose:
+ print('Numerical integration of ODE using implicite first-order method (Euler) was successful.')
+ return x, iterations
+
+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 t=0.
+ xp[0,:] = xp0 # Initial condition gives solution at t=0.
+ xpp[0,:] = xpp0 # Initial condition gives solution at t=0.
+ 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))
+ #x11 = x[i,:] + dxdt*Dt # Approximate solution at next value of x
+
+ 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)
+
+def newmark_newtonraphson_rdk(fnm, x0, xp0, xpp0, t, *p, gamma=.5, beta=.25, maxIterations=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 t=0.
+ xp[0,:] = xp0 # Initial condition gives solution at t=0.
+ xpp[0,:] = xpp0 # Initial condition gives solution at t=0.
+ 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(maxIterations): # fixed-point iteration
+ #dxdt = array(f(t[i+1], x1, p))
+ #x11 = x[i,:] + dxdt*Dt # Approximate solution at next value of x
+
+ 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)
diff --git a/src/solver_model.py b/src/solver_model.py
new file mode 100644
index 0000000..8e76a6c
--- /dev/null
+++ b/src/solver_model.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Mathmatical models governed by ordinary differential equations.
+
+Describes initial value problems as systems of first order ordinary differential
+equations.
+
+:Date: 2019-05-25
+
+.. module:: solver_model
+ :platform: *nix, Windows
+ :synopsis: Mechanical models.
+
+.. moduleauthor:: Daniel Weschke
+"""
+from __future__ import division, print_function
+from numpy import array, cos, sin, dot, square
+from numpy.linalg import inv
+
+def disk(x, t, *p):
+ """Rotation of an eccentric disk.
+
+ :param x: values of the function
+ :type x: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function
+
+ * diameter
+ * eccentricity
+ * torque
+ """
+ qp1 = x[3]
+ qp2 = x[4]
+ qp3 = x[5]
+ M = array([[1, 0, cos(x[2])], [0, 1, -sin(x[2])], [0, 0, 1]])
+ y = array([[-2*p[0]*x[3]+sin(x[2])*x[5]**2-2*p[0]*cos(x[2])*x[5]-x[0]], \
+ [-2*p[0]*x[4]+cos(x[2])*x[5]**2-2*p[0]*sin(x[2])*x[5]-x[1]], \
+ [p[2]-p[1]*x[1]*sin(x[2])+p[1]*x[0]*cos(x[2])]])
+ qp46 = dot(inv(M), y)
+ qp4, qp5, qp6 = qp46.reshape(-1,).tolist() # 2d array to 1d array to list
+ return qp1, qp2, qp3, qp4, qp5, qp6
+
+def disk_nm(xn, xpn, xppn, t, *p):
+ """Rotation of an eccentric disk.
+
+ :param xn: values of the function
+ :type xn: list
+ :param xpn: first derivative values of the function
+ :type xpn: list
+ :param xppn: second derivative values of the function
+ :type xppn: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function
+
+ * diameter
+ * eccentricity
+ * torque
+ """
+ N = array([[xppn[0]+cos(xn[2])*xppn[2]+2*p[0]*xpn[0]+2*p[0]*cos(xn[2])*xpn[2]-sin(xn[2])*square(xpn[2])+xn[0]],
+ [xppn[1]-sin(xn[2])*xppn[2]+2*p[0]*xpn[1]-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])+xn[1]],
+ [xppn[2]+p[1]*(-cos(xn[2])*xn[0]+sin(xn[2])*xn[1])-p[2]]])
+ dN = array([[1, 0, -sin(xn[2]*xppn[2])-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])],
+ [0, 1, -cos(xn[2]*xppn[2])-2*p[0]*cos(xn[2])*xpn[2]+sin(xn[2])*square(xpn[2])],
+ [-p[1]*cos(xn[2]), p[1]*cos(xn[2]), p[1]*(sin(xn[2])*xn[0]+cos(xn[2])*xn[1])]])
+ dNp = array([[2*p[0], 0, 2*p[0]*cos(xn[2])-2*sin(xn[2])*xpn[2]],
+ [0, 2*p[0], -2*p[0]*sin(xn[2])-2*cos(xn[2])*xpn[2]],
+ [0, 0, 0]])
+ dNpp = array([[1, 0, cos(xn[2])],
+ [0, 1, -sin(xn[2])],
+ [0, 0, 1]])
+ return N, dN, dNp, dNpp
+
+def disk_nmmdk(xn, xpn, xppn, t, *p):
+ """Rotation of an eccentric disk.
+
+ :param xn: values of the function
+ :type xn: list
+ :param xpn: derivative values of the function
+ :type xpn: list
+ :param xppn: second derivative values of the function
+ :type xppn: list
+ :param t: time
+ :type t: list
+ :param `*p`: parameters of the function
+
+ * diameter
+ * eccentricity
+ * torque
+ """
+ rm = array([[xppn[0]+cos(xn[2])*xppn[2]],
+ [xppn[1]-sin(xn[2])*xppn[2]],
+ [xppn[2]]])
+ rmx = array([[0, 0, -sin(xn[2]*xppn[2])],
+ [0, 0, -cos(xn[2]*xppn[2])],
+ [0, 0, 0]])
+ rmxpp = array([[1, 0, cos(xn[2])],
+ [0, 1, -sin(xn[2])],
+ [0, 0, 1]])
+ rd = array([[2*p[0]*xpn[0]+2*p[0]*cos(xn[2])*xpn[2]-sin(xn[2])*square(xpn[2])],
+ [2*p[0]*xpn[1]-2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])],
+ [0]])
+ rdx = array([[0, 0, -2*p[0]*sin(xn[2])*xpn[2]-cos(xn[2])*square(xpn[2])],
+ [0, 0, -2*p[0]*cos(xn[2])*xpn[2]+sin(xn[2])*square(xpn[2])],
+ [0, 0, 0]])
+ rdxp = array([[2*p[0], 0, 2*p[0]*cos(xn[2])-2*sin(xn[2])*xpn[2]],
+ [0, 2*p[0], -2*p[0]*sin(xn[2])-2*cos(xn[2])*xpn[2]],
+ [0, 0, 0]])
+ rk = array([[xn[0]],
+ [xn[1]],
+ [p[1]*(-cos(xn[2])*xn[0]+sin(xn[2])*xn[1])]])
+ rkx = array([[ 1, 0, 0],
+ [ 0, 1, 0],
+ [-p[1]*cos(xn[2]), p[1]*cos(xn[2]), p[1]*(sin(xn[2])*xn[0]+cos(xn[2])*xn[1])]])
+ f = array([[0], [0], [p[2]]])
+ return rm, rmx, rmxpp, rd, rdx, rdxp, rk, rkx, f
+
+if __name__ == '__main__':
+ 0
diff --git a/tests/test_fit.py b/tests/test_fit.py
new file mode 100644
index 0000000..fa7a7d0
--- /dev/null
+++ b/tests/test_fit.py
@@ -0,0 +1,16 @@
+import unittest
+
+import os
+import sys
+sys.path.insert(0, os.path.abspath('../src'))
+import fit
+
+
+class TestFit(unittest.TestCase):
+
+ def test_property(self):
+ self.assertEqual()
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/tests/test_solver.py b/tests/test_solver.py
new file mode 100644
index 0000000..fe41a83
--- /dev/null
+++ b/tests/test_solver.py
@@ -0,0 +1,201 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""Numerical solver.
+
+:Date: 2019-05-25
+
+.. module:: solver
+ :platform: *nix, Windows
+ :synopsis: Numerical solver.
+
+.. moduleauthor:: Daniel Weschke
+"""
+from __future__ import division, print_function
+import unittest
+import math
+from numpy import linspace, allclose
+from numpy.linalg import inv
+
+import os
+import sys
+sys.path.insert(0, os.path.abspath('../src'))
+import solver
+import solver_model
+
+#### Post processing
+from matplotlib.pyplot import figure, subplots, plot, show
+from numpy import sqrt, square, concatenate
+
+def plotx1rphi(x, t, title):
+ """Plot plane rotation (x, y, phi)
+ """
+ x1 = x[:, 0]
+ x2 = x[:, 1]
+ x3 = x[:, 2]
+ xp1 = x[:, 3]
+ xp2 = x[:, 4]
+ xp3 = x[:, 5]
+ fig, ax1 = subplots()
+ ax1.set_title(title)
+ ax1.plot(t, x1, 'b-')
+ #ax1.plot(t, x3, 'b-')
+ ax1.set_xlabel('time t in s')
+ # Make the y-axis label and tick labels match the line color.
+ ax1.set_ylabel('x_1', color='b')
+ for tl in ax1.get_yticklabels():
+ tl.set_color('b')
+ envelope = sqrt(square(x1)+square(x2))
+ ax1.plot(t, envelope, 'b--', t, -envelope, 'b--')
+ ax2 = ax1.twinx()
+ ax2.plot(t, xp3, 'r')
+ ax2.set_ylabel('phi', color='r')
+ for tl in ax2.get_yticklabels():
+ tl.set_color('r')
+ if max(xp3) < 2:
+ ax2.set_ylim([0,2])
+
+f = solver_model.disk
+fnm = solver_model.disk_nm
+fnmmdk = solver_model.disk_nmmdk
+
+
+class TestSolver(unittest.TestCase):
+
+ Dt = .01 # Delta t
+ tend = .1 # t_end
+ #tend = 300
+ nt = math.ceil(tend/Dt) # number of timesteps
+ t = linspace(0., tend, nt+1) # time vector
+ x0 = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 # initial conditions
+ p = .01, .002, .015 # parameter of the model
+
+ def test_e1(self):
+ xe1 = solver.e1(f, self.x0, self.t, *self.p)
+ #plotx1rphi(xe1, self.t, ("Euler (Runge-Kutta 1th order), increments: %d" % self.nt))
+ #show()
+ self.assertTrue(allclose(xe1[:, 5],
+ [0., 0.00015, 0.0003, 0.00045, 0.0006, 0.00075, 0.0009, 0.00105, 0.0012,
+ 0.00135, 0.0015 ]))
+
+ Dt_2 = Dt/5.
+ nt_2 = math.ceil(tend/Dt_2)
+ t_2 = linspace(0., tend, nt_2+1)
+
+ def test_e1_2(self):
+ xe1_2 = solver.e1(f, self.x0, self.t_2, *self.p)
+ #plotx1rphi(xe1_2, self.t_2, ("Euler (Runge-Kutta 1th order), increments: %d" % self.nt_2))
+ #show()
+ self.assertTrue(allclose(xe1_2[:, 5],
+[0.00000000e+00, 3.00000000e-05, 6.00000000e-05, 8.99999998e-05,
+ 1.19999999e-04, 1.49999998e-04, 1.79999995e-04, 2.09999992e-04,
+ 2.39999987e-04, 2.69999980e-04, 2.99999971e-04, 3.29999960e-04,
+ 3.59999947e-04, 3.89999931e-04, 4.19999913e-04, 4.49999891e-04,
+ 4.79999866e-04, 5.09999837e-04, 5.39999804e-04, 5.69999767e-04,
+ 5.99999726e-04, 6.29999681e-04, 6.59999630e-04, 6.89999575e-04,
+ 7.19999514e-04, 7.49999448e-04, 7.79999376e-04, 8.09999298e-04,
+ 8.39999214e-04, 8.69999123e-04, 8.99999026e-04, 9.29998921e-04,
+ 9.59998810e-04, 9.89998691e-04, 1.01999856e-03, 1.04999843e-03,
+ 1.07999829e-03, 1.10999814e-03, 1.13999798e-03, 1.16999781e-03,
+ 1.19999763e-03, 1.22999744e-03, 1.25999725e-03, 1.28999704e-03,
+ 1.31999682e-03, 1.34999660e-03, 1.37999636e-03, 1.40999611e-03,
+ 1.43999585e-03, 1.46999558e-03, 1.49999530e-03]
+ ))
+
+ maxIterations = 1000
+ tol = 1e-9 # against residuum
+
+ def test_e2(self):
+ xe2 = solver.e2(f, self.x0, self.t, *self.p)
+ #plotx1rphi(xe2, self.t, ("Runge-Kutta 2th order, increments: %d" % self.nt))
+ #show()
+ self.assertTrue(allclose(xe2[:, 5],
+ [0., 0.00015, 0.0003, 0.00045, 0.0006, 0.00075, 0.0009, 0.00105, 0.0012,
+ 0.00135, 0.0015 ]))
+
+ def test_e4(self):
+ xe4 = solver.e4(f, self.x0, self.t, *self.p)
+ #plotx1rphi(xe4, self.t, ("Runge-Kutta 4th order, increments: %d" % self.nt))
+ #show()
+ self.assertTrue(allclose(xe4[:, 5],
+ [0., 0.00015, 0.0003, 0.00045, 0.0006, 0.00075, 0.0009, 0.00105, 0.0012,
+ 0.00135, 0.0015 ]))
+
+ def test_i1(self):
+ xi1, iterations = solver.i1(f, self.x0, self.t, *self.p, self.maxIterations, self.tol)
+ #plotx1rphi(xi1, self.t, ("Backward Euler, increments: %d, iterations: %d" % (self.nt, sum(iterations))))
+ #show()
+ self.assertTrue(allclose(xi1[:, 5],
+ [0., 0.00015, 0.0003, 0.00045, 0.0006, 0.00075, 0.0009, 0.00105, 0.0012,
+ 0.00135, 0.0015 ]))
+
+ gamma = .5 # newmark parameter for velocity
+ beta = .25 # newmark parameter for displacement
+
+ def test_nm(self):
+ xnm, xpnm, xppnm, iterations = solver.newmark_newtonraphson(fnm,
+ (0,0,0), (0,0,0), (0,0,0), self.t, *self.p,
+ self.gamma, self.beta, 1, self.tol)
+ #print(xnm)
+ #print(xpnm)
+ x = concatenate((xnm, xpnm), axis=1)
+ #plotx1rphi(x, self.t, ("Newmark, increments: %d, iterations: %d" % (self.nt, sum(iterations))))
+ #show()
+ self.assertTrue(allclose(x[:, 5],
+[0.00000000e+00, 7.49999925e-05, 2.24999951e-04, 3.74999839e-04,
+ 5.24999621e-04, 6.74999269e-04, 8.24998752e-04, 9.74998039e-04,
+ 1.12499710e-03, 1.27499591e-03, 1.42499444e-03]
+ ))
+
+ def test_nm_2(self):
+ xnm, xpnm, xppnm, iterations = solver.newmark_newtonraphson(fnm,
+ (0,0,0), (0,0,0), (0,0,0), self.t, *self.p,
+ self.gamma, self.beta, self.maxIterations, self.tol)
+ #print(xnm)
+ #print(xpnm)
+ x = concatenate((xnm, xpnm), axis=1)
+ #plotx1rphi(x, self.t, ("Newmark, increments: %d, iterations: %d" % (self.nt, sum(iterations))))
+ #show()
+ self.assertTrue(allclose(x[:, 5],
+[0.00000000e+00, 7.49999925e-05, 2.24999951e-04, 3.74999839e-04,
+ 5.24999621e-04, 6.74999269e-04, 8.24998752e-04, 9.74998039e-04,
+ 1.12499710e-03, 1.27499591e-03, 1.42499444e-03]
+ ))
+
+ def test_nmmdk(self):
+ xnm, xpnm, xppnm, iterations = solver.newmark_newtonraphson_rdk(fnmmdk,
+ (0,0,0), (0,0,0), (0,0,0), self.t, *self.p,
+ self.gamma, self.beta, 1, self.tol)
+ #print(xnm)
+ #print(xpnm)
+ x = concatenate((xnm, xpnm), axis=1)
+ #plotx1rphi(x, self.t, ("Newmark, increments: %d, iterations: %d" % (self.nt, sum(iterations))))
+ #show()
+ self.assertTrue(allclose(x[:, 5],
+[0.00000000e+00, 7.49999963e-05, 2.24999974e-04, 3.74999906e-04,
+ 5.24999764e-04, 6.74999516e-04, 8.24999134e-04, 9.74998586e-04,
+ 1.12499784e-03, 1.27499688e-03, 1.42499566e-03]
+ ))
+
+ def test_nmmdk_2(self):
+ xnm, xpnm, xppnm, iterations = solver.newmark_newtonraphson_rdk(fnmmdk,
+ (0,0,0), (0,0,0), (0,0,0), self.t, *self.p,
+ self.gamma, self.beta, 100, self.tol)
+ #print(xnm)
+ #print(xpnm)
+ x = concatenate((xnm, xpnm), axis=1)
+ #plotx1rphi(x, self.t, ("Newmark, increments: %d, iterations: %d" % (self.nt, sum(iterations))))
+ #show()
+ self.assertTrue(allclose(x[:, 5],
+[0.00000000e+00, 7.49999963e-05, 2.24999974e-04, 3.74999906e-04,
+ 5.24999764e-04, 6.74999516e-04, 8.24999134e-04, 9.74998586e-04,
+ 1.12499784e-03, 1.27499688e-03, 1.42499566e-03]
+ ))
+
+ #from scipy.integrate import odeint
+ #xode = odeint(solver_model.disk, x0, t, p, printmessg=True)
+ #plotx1rphi(xode, t, ("ODE, Inkremente: %d" % nt))
+ #show()
+
+
+if __name__ == '__main__':
+ unittest.main()