From d0873a36da06ca325d47a29644f0007a449bb7e7 Mon Sep 17 00:00:00 2001 From: Daniel Weschke Date: Sun, 26 May 2019 21:54:58 +0200 Subject: [PATCH] add ode solver, tests and update docs --- docs/build/html/.buildinfo | 2 +- docs/build/html/_sources/modules.rst.txt | 2 + docs/build/html/_sources/solver.rst.txt | 7 + docs/build/html/_sources/solver_model.rst.txt | 7 + docs/build/html/_sources/src.rst.txt | 46 --- docs/build/html/_static/alabaster.css | 4 +- docs/build/html/data.html | 50 ++- docs/build/html/date.html | 99 +++-- docs/build/html/fit.html | 47 ++- docs/build/html/genindex.html | 55 ++- docs/build/html/geometry.html | 188 +++++++-- docs/build/html/index.html | 1 + docs/build/html/modules.html | 3 + docs/build/html/objects.inv | 13 +- docs/build/html/py-modindex.html | 33 +- docs/build/html/search.html | 1 + docs/build/html/searchindex.js | 2 +- docs/build/html/solver.html | 289 ++++++++++++++ docs/build/html/solver_model.html | 177 +++++++++ docs/build/html/src.html | 119 ------ docs/source/conf.py | 2 + docs/source/modules.rst | 2 + docs/source/solver.rst | 7 + docs/source/solver_model.rst | 7 + docs/source/src.rst | 46 --- src/data.py | 43 ++- src/date.py | 82 ++-- src/fit.py | 48 ++- src/geometry.py | 141 +++++-- src/solver.py | 357 ++++++++++++++++++ src/solver_model.py | 120 ++++++ tests/test_fit.py | 16 + tests/test_solver.py | 201 ++++++++++ 33 files changed, 1821 insertions(+), 396 deletions(-) create mode 100644 docs/build/html/_sources/solver.rst.txt create mode 100644 docs/build/html/_sources/solver_model.rst.txt delete mode 100644 docs/build/html/_sources/src.rst.txt create mode 100644 docs/build/html/solver.html create mode 100644 docs/build/html/solver_model.html delete mode 100644 docs/build/html/src.html create mode 100644 docs/source/solver.rst create mode 100644 docs/source/solver_model.rst delete mode 100644 docs/source/src.rst create mode 100644 src/solver.py create mode 100644 src/solver_model.py create mode 100644 tests/test_fit.py create mode 100644 tests/test_solver.py diff --git a/docs/build/html/.buildinfo b/docs/build/html/.buildinfo index 8c381aa..8106eaa 100644 --- a/docs/build/html/.buildinfo +++ b/docs/build/html/.buildinfo @@ -1,4 +1,4 @@ # Sphinx build info version 1 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. -config: 09cf9ce4470e7ee74b6ce2abe5e4453a +config: 94375f4299332632f508e2242b7c30d8 tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/docs/build/html/_sources/modules.rst.txt b/docs/build/html/_sources/modules.rst.txt index 8b02d0a..6a908e5 100644 --- a/docs/build/html/_sources/modules.rst.txt +++ b/docs/build/html/_sources/modules.rst.txt @@ -8,3 +8,5 @@ src date fit geometry + solver + solver_model diff --git a/docs/build/html/_sources/solver.rst.txt b/docs/build/html/_sources/solver.rst.txt new file mode 100644 index 0000000..0d6d757 --- /dev/null +++ b/docs/build/html/_sources/solver.rst.txt @@ -0,0 +1,7 @@ +solver module +============= + +.. automodule:: solver + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/build/html/_sources/solver_model.rst.txt b/docs/build/html/_sources/solver_model.rst.txt new file mode 100644 index 0000000..bc08685 --- /dev/null +++ b/docs/build/html/_sources/solver_model.rst.txt @@ -0,0 +1,7 @@ +solver\_model module +==================== + +.. automodule:: solver_model + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/build/html/_sources/src.rst.txt b/docs/build/html/_sources/src.rst.txt deleted file mode 100644 index aa0f83a..0000000 --- a/docs/build/html/_sources/src.rst.txt +++ /dev/null @@ -1,46 +0,0 @@ -src package -=========== - -Submodules ----------- - -src.data module ---------------- - -.. automodule:: src.data - :members: - :undoc-members: - :show-inheritance: - -src.date module ---------------- - -.. automodule:: src.date - :members: - :undoc-members: - :show-inheritance: - -src.fit module --------------- - -.. automodule:: src.fit - :members: - :undoc-members: - :show-inheritance: - -src.geometry module -------------------- - -.. automodule:: src.geometry - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: src - :members: - :undoc-members: - :show-inheritance: diff --git a/docs/build/html/_static/alabaster.css b/docs/build/html/_static/alabaster.css index e127322..b1d8d1a 100644 --- a/docs/build/html/_static/alabaster.css +++ b/docs/build/html/_static/alabaster.css @@ -310,8 +310,8 @@ div.hint { } div.seealso { - background-color: #EEE; - border: 1px solid #CCC; + background-color: #3c3c3c; + border: 1px solid #2C2C2C; } div.topic { diff --git a/docs/build/html/data.html b/docs/build/html/data.html index 7525a6b..bf20189 100644 --- a/docs/build/html/data.html +++ b/docs/build/html/data.html @@ -13,6 +13,7 @@ + @@ -33,28 +34,67 @@

data module

-
+

Read and write data to or from file.

+
data_load(file_name, verbose=False)
-

load stored program objects from binary file

+

Load stored program objects from binary file.

+
+
Parameters
+
    +

  • file_name (str) – file to load

    • +

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

    • +
+
+
Returns
+

loaded data

+
+
Return type
+

object

+
+
data_read(file_name, x_column, y_column)
-

read ascii data file

+

Read ascii data file.

+
+
Parameters
+
    +

  • filename (str) – file to read

    • +

  • x_column (int) – column index for the x data (first column is 0)

    • +

  • y_column (int) – column index for the y data (first column is 0)

    • +
+
+
Returns
+

x and y

+
+
Return type
+

tuple(list, list)

+
+
data_store(file_name, object_data)
-

store program objects to binary file

+

Store program objects to binary file.

+
+
Parameters
+
    +

  • file_name (str) – file to store

    • +

  • object_data (object) – data to store

    • +
+
+
main()
-
+

Main function.

+
diff --git a/docs/build/html/date.html b/docs/build/html/date.html index 1edfb8e..e0fb1b6 100644 --- a/docs/build/html/date.html +++ b/docs/build/html/date.html @@ -13,6 +13,7 @@ + @@ -33,15 +34,25 @@

date module

-

Created on Mon Jan 15 21:52:34 2018

-

@author: Daniel Weschke

-
+

Calculate spacial dates.

+
+
Date
+

2018-01-15

+
+
+
ascension_of_jesus(year)

Ascension of Jesus.

-
Returns
-

datetime.date – the day of ascension of Jesus

+
Parameters
+

year (int) – the year to calculate the ascension of Jesus

+
+
Returns
+

the day of ascension of Jesus

+
+
Return type
+

datetime.date

@@ -51,8 +62,14 @@ easter_friday(year)

Easter Friday.

-
Returns
-

datetime.date – the day of Easter Friday

+
Parameters
+

year (int) – the year to calculate the Easter Friday

+
+
Returns
+

the day of Easter Friday

+
+
Return type
+

datetime.date

@@ -62,8 +79,14 @@ easter_monday(year)

Easter Monday.

-
Returns
-

datetime.date – the day of Easter Monday

+
Parameters
+

year (int) – the year to calculate the Easter Monday

+
+
Returns
+

the day of Easter Monday

+
+
Return type
+

datetime.date

@@ -73,8 +96,14 @@ easter_sunday(year)

Easter Sunday.

-
Returns
-

datetime.date – the day of Easter Sunday

+
Parameters
+

year (int) – the year to calculate the Easter Sunday

+
+
Returns
+

the day of Easter Sunday

+
+
Return type
+

datetime.date

@@ -84,24 +113,32 @@ gaußsche_osterformel(year)

Gaußsche Osterformel.

-
Returns
-

int – the day of Easter Sunday as a day in march.

+
Parameters
+

year (int) – the year to calculate the Easter Sunday

+
+
Returns
+

the day of Easter Sunday as a day in march.

+
+
Return type
+

int

+
+
Variables
+
    +

  • X (int) – Das Jahr / year

    • +

  • K(X) (int) – Die Säkularzahl

    • +

  • M(X) (int) – Die säkulare Mondschaltung

    • +

  • S(K) (int) – Die säkulare Sonnenschaltung

    • +

  • A(X) (int) – Den Mondparameter

    • +

  • D(A,M) (int) – Den Keim für den ersten Vollmond im Frühling

    • +

  • R(D,A) (int) – Die kalendarische Korrekturgröße

    • +

  • OG(D,R) (int) – Die Ostergrenze

    • +

  • SZ(X,S) (int) – Den ersten Sonntag im März

    • +

  • OE(OG,SZ) (int) – Die Entfernung des Ostersonntags von der Ostergrenze (Osterentfernung in Tagen)

    • +

  • OS(OG,OE) (int) – Das Datum des Ostersonntags als Märzdatum (32. März = 1. April usw.)

    • +

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

@@ -110,8 +147,14 @@ pentecost(year)

Pentecost.

-
Returns
-

datetime.date – the day of Pentecost

+
Parameters
+

year (int) – the year to calculate the Pentecost

+
+
Returns
+

the day of Pentecost

+
+
Return type
+

datetime.date

diff --git a/docs/build/html/fit.html b/docs/build/html/fit.html index e1f1896..e0b6e50 100644 --- a/docs/build/html/fit.html +++ b/docs/build/html/fit.html @@ -13,6 +13,7 @@ + @@ -33,26 +34,32 @@

fit module

-
+

Function and approximation.

+
gauss(x, *p)

Gauss distribution function.

+
+\[f(x)=ae^{-(x-b)^{2}/(2c^{2})}\]
Parameters
    -

  • x – positions where the gauss function will be calculated

    • -

  • p –

    gauss parameters [a, b, c, d]:

    +

  • x (int or float or list or numpy.ndarray) – positions where the gauss function will be calculated

    • +

  • p (list) –

    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 (\(\int y \,\mathrm{d}x=1 \Leftrightarrow a=1/(c\sqrt{2\pi})\) )

      • +

    • b – expected value \(\mu\) (position of maximum, default = 0)

      • +

    • c – standard deviation \(\sigma\) (variance \(\sigma^2=c^2\))

    • d – vertical offset (default = 0)

Returns
-

array – gauss values at given positions x

+

gauss values at given positions x

+
+
Return type
+

numpy.ndarray

@@ -64,22 +71,32 @@
Parameters
    -

  • x – positions

    • -

  • y – values

    • -

  • e – error values

    • -

  • x_fit – positions of fitted function (default steps is 3*len(x) but min 150)

    • +

  • x (int or float or list or numpy.ndarray) – positions

    • +

  • y (int or float or list or numpy.ndarray) – values

    • +

  • e (int or float or list or numpy.ndarray) – error values (default = None)

    • +

  • x_fit (int or float or list or numpy.ndarray) – positions of fitted function (default = None, if None then x +is used)

    • +

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

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 \(\mu\), standard deviation +\(\sigma\), vertical offset d)

    • +

  • numpy.float64 – full width at half maximum (FWHM)

+
Return type
+

tuple

+
+
+

See also

+

gauss()

+
diff --git a/docs/build/html/genindex.html b/docs/build/html/genindex.html index ebd1f99..423ac4b 100644 --- a/docs/build/html/genindex.html +++ b/docs/build/html/genindex.html @@ -14,6 +14,7 @@ + @@ -42,8 +43,10 @@ | E | F | G + | I | L | M + | N | P | R | S @@ -73,17 +76,23 @@

D

- +
@@ -91,10 +100,16 @@

E

+
+ +

I

+ +
@@ -146,6 +169,18 @@ +

N

+ + + +
+

P

    @@ -175,6 +210,12 @@

    S

    + diff --git a/docs/build/html/geometry.html b/docs/build/html/geometry.html index 50bca07..0f3c4fa 100644 --- a/docs/build/html/geometry.html +++ b/docs/build/html/geometry.html @@ -13,6 +13,7 @@ + @@ -33,32 +34,81 @@

    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

    +
    +
    +
    +

    See also

    +

    cubic()

    +
    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
    +
      +

    • width (int or float) – the width of the rectangle

      • +

    • height (int or float) – the height of the rectangle

      • +
    +
    +
    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
    +
      +

    • origin (tuple) – the center of rotation

      • +

    • angle (int or float) – the rotation angle

      • +

    • *pts – points to rotate

      • +

    • **kwargs – options

      • +
    +
    +
    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
    +
      +

    • origin (tuple) – the center of rotation

      • +

    • angle (int or float) – the rotation angle

      • +

    • *pts – points to rotate

      • +

    • **kwargs – options

      • +
    +
    +
    Returns
    +

    (point_x, point_y) or (point1, point2, …)

    +
    +
    Return type
    +

    tuple

    +
    +

    See also

    +

    rotate()

    +
    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

    +
    +
    +
    +

    See also

    +

    rectangle()

    +
    translate(vec, *pts)

    Translate a point or polygon by a given vector.

    -
    -
    returns

    (point_x, point_y) -or -(point1, point2, …)

    +
    +
    Parameters
    +
      +

    • vec (tuple) – translation vector

      • +

    • *pts – points to translate

      • +
    +
    +
    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 @@ + diff --git a/docs/build/html/modules.html b/docs/build/html/modules.html index d3754f5..1d347e6 100644 --- a/docs/build/html/modules.html +++ b/docs/build/html/modules.html @@ -13,6 +13,7 @@ + @@ -39,6 +40,8 @@
  • 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!ZtcvJ(aGa т]QYФ -Cȉ0!?XqS -G -*< 4 -?<؊ш?Nn> -D4swD71QݹɼoMQ)NCtKӴj0YT5"ϟ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Ӡz0xu ϋ2XX wa +[0X ޝM.|8H>-UK +%J{+ztjC +g*lm m3jV}Wr +-x 3sPʼ~hh8њ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 @@ + @@ -44,7 +45,8 @@
    d | f | - g + g | + s
    @@ -54,29 +56,42 @@ + data(*nix, Windows) + date(*nix, Windows) + fit(*nix, Windows) + geometry(*nix, Windows) + + + + + + + +
    - data -
    + Handle data files.
    - date -
    + Special dates.
     
    f
    - fit -
    + Function and approximation.
     
    g
    - geometry -
    + 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 @@ + diff --git a/docs/build/html/searchindex.js b/docs/build/html/searchindex.js index 6da2a79..5c989c8 100644 --- a/docs/build/html/searchindex.js +++ b/docs/build/html/searchindex.js @@ -1 +1 @@ -Search.setIndex({docnames:["data","date","fit","geometry","index","modules","src"],envversion:{"sphinx.domains.c":1,"sphinx.domains.changeset":1,"sphinx.domains.cpp":1,"sphinx.domains.javascript":1,"sphinx.domains.math":2,"sphinx.domains.python":1,"sphinx.domains.rst":1,"sphinx.domains.std":1,sphinx:56},filenames:["data.rst","date.rst","fit.rst","geometry.rst","index.rst","modules.rst","src.rst"],objects:{"":{data:[0,0,0,"-"],date:[1,0,0,"-"],fit:[2,0,0,"-"],geometry:[3,0,0,"-"]},data:{data_load:[0,1,1,""],data_read:[0,1,1,""],data_store:[0,1,1,""],main:[0,1,1,""]},date:{"gau\u00dfsche_osterformel":[1,1,1,""],ascension_of_jesus:[1,1,1,""],easter_friday:[1,1,1,""],easter_monday:[1,1,1,""],easter_sunday:[1,1,1,""],pentecost:[1,1,1,""]},fit:{gauss:[2,1,1,""],gauss_fit:[2,1,1,""],main:[2,1,1,""]},geometry:{cubic:[3,1,1,""],cubic_deg:[3,1,1,""],line:[3,1,1,""],points:[3,1,1,""],rectangle:[3,1,1,""],rotate:[3,1,1,""],rotate_deg:[3,1,1,""],square:[3,1,1,""],translate:[3,1,1,""]}},objnames:{"0":["py","module","Python module"],"1":["py","function","Python function"]},objtypes:{"0":"py:module","1":"py:function"},terms:{"9fsche_osterformel":1,"default":2,"f\u00fcr":1,"fr\u00fchling":1,"function":2,"gau\u00dfsch":1,"gau\u00dfsche_osterformel":1,"int":1,"korrekturgr\u00f6\u00df":1,"m\u00e4rz":1,"m\u00e4rzdatum":1,"return":[1,2,3],"s\u00e4kular":1,"s\u00e4kularzahl":1,Das:1,The:3,algorithmu:1,als:1,amplitud:2,angl:3,angle1:3,angle2:3,april:1,around:3,arrai:2,ascens:1,ascension_of_jesu:1,ascii:0,author:[1,3],binari:0,calcul:2,counterclockwis:3,creat:1,cubic:3,cubic_deg:3,dai:1,daniel:[1,3],data:[2,5],data_load:0,data_read:0,data_stor:0,date:[3,5],datetim:1,datum:1,degre:3,den:1,der:1,des:1,deviat:2,die:1,distribut:2,easter:1,easter_fridai:1,easter_mondai:1,easter_sundai:1,entfernung:1,error:2,ersten:1,expect:2,fals:[0,2],file:0,file_nam:0,fit:5,frac:[],fridai:1,from:0,full:2,fwhm:2,gau:1,gauss:2,gauss_fit:2,geometri:5,gilt:1,given:[2,3],gregorianischen:1,half:2,height:3,http:1,index:4,integr:2,jahr:1,jan:1,jesu:1,kalend:1,kalendarisch:1,keim:1,kwarg:3,len:2,line:3,load:0,main:[0,2],march:1,maximum:2,min:2,modul:[4,5],mon:1,mondai:1,mondparamet:1,mondschaltung:1,none:2,number:3,object:0,object_data:0,offset:2,org:1,origin:3,osterentfernung:1,osterformel:1,ostergrenz:1,ostersonntag:1,page:4,paramet:2,pentecost:1,point1:3,point1_i:3,point1_x:3,point2:3,point2_i:3,point2_x:3,point3:3,point4:3,point:3,point_i:3,point_x:3,points1_i:3,polygon:3,popt:2,posit:2,program:0,pts:3,rac:[],radian:3,read:0,rectangl:3,rotat:3,rotate_deg:3,sampl:3,sample_point1_x:3,sample_point2_i:3,sample_point2_x:3,sample_points1_i:3,search:4,should:3,sigma:2,sonnenschaltung:1,sonntag:1,sourc:1,sqrt:2,squar:3,standard:2,step:2,store:0,sum_:[],sundai:1,tagen:1,test:2,translat:3,usw:1,valu:2,varianc:2,vec:3,vector:3,verbos:[0,2],vertic:2,vollmond:1,von:1,weschk:[1,3],where:2,width:[2,3],wiki:1,wikipedia:1,x_column:0,x_fit:2,y_column:0,y_fit:2,year:1},titles:["data module","date module","fit module","geometry module","Welcome to pylib\u2019s documentation!","src","src package"],titleterms:{content:6,data:[0,6],date:[1,6],document:4,fit:[2,6],geometri:[3,6],indic:4,modul:[0,1,2,3,6],packag:6,pylib:4,src:[5,6],submodul:6,tabl:4,welcom:4}}) \ No newline at end of file +Search.setIndex({docnames:["data","date","fit","geometry","index","modules","solver","solver_model"],envversion:{"sphinx.domains.c":1,"sphinx.domains.changeset":1,"sphinx.domains.cpp":1,"sphinx.domains.javascript":1,"sphinx.domains.math":2,"sphinx.domains.python":1,"sphinx.domains.rst":1,"sphinx.domains.std":1,sphinx:56},filenames:["data.rst","date.rst","fit.rst","geometry.rst","index.rst","modules.rst","solver.rst","solver_model.rst"],objects:{"":{data:[0,0,0,"-"],date:[1,0,0,"-"],fit:[2,0,0,"-"],geometry:[3,0,0,"-"],solver:[6,0,0,"-"],solver_model:[7,0,0,"-"]},data:{data_load:[0,1,1,""],data_read:[0,1,1,""],data_store:[0,1,1,""],main:[0,1,1,""]},date:{"gau\u00dfsche_osterformel":[1,1,1,""],ascension_of_jesus:[1,1,1,""],easter_friday:[1,1,1,""],easter_monday:[1,1,1,""],easter_sunday:[1,1,1,""],pentecost:[1,1,1,""]},fit:{gauss:[2,1,1,""],gauss_fit:[2,1,1,""],main:[2,1,1,""]},geometry:{cubic:[3,1,1,""],cubic_deg:[3,1,1,""],line:[3,1,1,""],points:[3,1,1,""],rectangle:[3,1,1,""],rotate:[3,1,1,""],rotate_deg:[3,1,1,""],square:[3,1,1,""],translate:[3,1,1,""]},solver:{e1:[6,1,1,""],e2:[6,1,1,""],e4:[6,1,1,""],i1:[6,1,1,""],newmark_newtonraphson:[6,1,1,""],newmark_newtonraphson_rdk:[6,1,1,""]},solver_model:{disk:[7,1,1,""],disk_nm:[7,1,1,""],disk_nmmdk:[7,1,1,""]}},objnames:{"0":["py","module","Python module"],"1":["py","function","Python function"]},objtypes:{"0":"py:module","1":"py:function"},terms:{"1st":6,"1th":[],"2nd":6,"4th":6,"9fsche_osterformel":1,"default":[0,2,6],"f\u00fcr":1,"float":[2,3,6],"fr\u00fchling":1,"function":[0,2,6,7],"gau\u00dfsch":1,"gau\u00dfsche_osterformel":1,"int":[0,1,2,3,6],"korrekturgr\u00f6\u00df":1,"m\u00e4rz":1,"m\u00e4rzdatum":1,"return":[0,1,2,3],"s\u00e4kular":1,"s\u00e4kularzahl":1,"vorw\u00e4rt":6,Das:1,One:6,The:[3,6],against:6,algorithmu:1,als:1,amplitud:2,angl:3,angle1:3,angle2:3,approxim:[2,6],april:1,around:3,ascens:1,ascension_of_jesu:1,ascii:0,backward:6,becom:6,begin:6,beta:6,binari:0,bmatrix:6,bool:[0,2,6],calcul:[1,2],cauchi:6,center:3,choos:6,column:0,condit:6,counterclockwis:3,cubic:3,cubic_deg:3,dai:1,data:[2,5],data_load:0,data_read:0,data_stor:0,date:[3,5,6,7],datetim:1,datum:1,ddot:6,degre:3,den:1,der:1,deriv:7,des:1,describ:7,deviat:2,diamet:[6,7],die:1,differenti:[6,7],disk:7,disk_nm:7,disk_nmmdk:7,displac:6,distribut:2,dot:6,easter:1,easter_fridai:1,easter_mondai:1,easter_sundai:1,eccentr:7,edg:3,end:[3,6],entfernung:1,equat:[6,7],error:[2,6],ersten:1,euler:6,everi:6,exampl:6,expect:2,explicit:6,explizit:6,fals:[0,2,6],file:0,file_nam:0,filenam:0,first:[0,6,7],fit:5,float64:2,fnm:6,forward:6,fourth:6,frac:3,fridai:1,from:[0,6],full:2,fwhm:2,gamma:6,gau:1,gauss:2,gauss_fit:2,geometri:5,gilt:1,given:[2,3,6],global:6,govern:7,gregorianischen:1,half:2,has:6,height:3,http:1,implicit:6,index:[0,4],inform:[0,2,6],initi:[6,7],iter:6,jahr:1,jesu:1,kalend:1,kalendarisch:1,keim:1,kutta:6,kwarg:3,ldot:6,leftrightarrow:2,line:3,list:[0,2,6,7],load:0,local:6,main:[0,2],march:1,mathmat:7,mathrm:2,max_iter:6,maximum:[2,6],maxiter:6,mean:6,mechan:[],method:6,model:7,modul:[4,5],mondai:1,mondparamet:1,mondschaltung:1,ndarrai:2,newmark:6,newmark_newtonraphson:6,newmark_newtonraphson_rdk:6,none:2,number:[3,6],numer:6,numpi:2,object:[0,3],object_data:0,ode1back:[],offset:2,one:[3,6],option:3,order:[6,7],ordinari:[6,7],org:1,origin:3,osterentfernung:1,osterformel:1,ostergrenz:1,ostersonntag:1,other:3,page:4,param:[],paramet:[0,1,2,3,6,7],pentecost:1,per:6,point1:3,point1_i:3,point1_x:3,point2:3,point2_i:3,point2_x:3,point3:3,point4:3,point:3,point_i:3,point_x:3,points1_i:3,polygon:3,polygonzugverfahren:6,popt:2,posit:2,print:6,printmessg:[],problem:[6,7],program:0,proport:6,pts:3,quad:6,radian:3,read:0,rearrang:3,rectangl:3,residuum:6,rk2:[],rk4:[],rotat:[3,7],rotate_deg:3,rung:6,sampl:3,sample_point1_x:3,sample_point2_i:3,sample_point2_x:3,sample_points1_i:3,sche:6,search:4,second:[6,7],set:6,should:3,sigma:2,size:[3,6],slope:3,solut:6,solv:6,solver:5,solver_model:5,sonnenschaltung:1,sonntag:1,sourc:1,spacial:1,sqrt:2,squar:[3,6],stabl:6,standard:[2,6],step:6,store:0,str:0,sundai:1,system:[6,7],t_0:6,t_i:6,tagen:1,test:2,thick:6,time:[6,7],tol:6,toler:6,torqu:7,translat:3,tupl:[0,2,3],type:[0,1,2,3],used:2,usw:1,valu:[2,6,7],variabl:1,varianc:2,vec:3,vector:3,veloc:6,verbos:[0,2,6],verfahren:6,vertic:2,vollmond:1,von:1,where:2,which:6,width:[2,3],wiki:1,wikipedia:1,write:0,x_0:6,x_1:[3,6],x_2:[3,6],x_column:0,x_fit:2,x_i:6,xp0:6,xpn:7,xpp0:6,xppn:7,y_1:3,y_2:3,y_column:0,y_fit:2,year:1},titles:["data module","date module","fit module","geometry module","Welcome to pylib\u2019s documentation!","src","solver module","solver_model module"],titleterms:{data:0,date:1,document:4,fit:2,geometri:3,indic:4,modul:[0,1,2,3,6,7],pylib:4,solver:6,solver_model:7,src:5,tabl:4,welcom:4}}) \ No newline at end of file diff --git a/docs/build/html/solver.html b/docs/build/html/solver.html new file mode 100644 index 0000000..8ef563e --- /dev/null +++ b/docs/build/html/solver.html @@ -0,0 +1,289 @@ + + + + + + + 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)

      • +
    +
    +
    +
    + +
    + + +
    + +
    +
    + +
    +
    +
    + ©2019, Daniel Weschke. + + | + Powered by Sphinx 2.0.1 + & Alabaster 0.7.12 + + | + Page source +
    + + + + + + \ No newline at end of file diff --git a/docs/build/html/solver_model.html b/docs/build/html/solver_model.html new file mode 100644 index 0000000..597c46c --- /dev/null +++ b/docs/build/html/solver_model.html @@ -0,0 +1,177 @@ + + + + + + + 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
    +
      +

    • x (list) – values of the function

      • +

    • t (list) – time

      • +

    • *p –

      parameters of the function

      +
        +

      • diameter

        • +

      • eccentricity

        • +

      • torque

        • +
      +

      • +
    +
    +
    +
    + +
    +
    +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

        • +
      +

      • +
    +
    +
    +
    + +
    + + +
    + +
    +
    + +
    +
    +
    + ©2019, Daniel Weschke. + + | + Powered by Sphinx 2.0.1 + & Alabaster 0.7.12 + + | + Page source +
    + + + + + + \ No newline at end of file diff --git a/docs/build/html/src.html b/docs/build/html/src.html deleted file mode 100644 index 1d9a351..0000000 --- a/docs/build/html/src.html +++ /dev/null @@ -1,119 +0,0 @@ - - - - - - - src package — pylib 2019.5.19 documentation - - - - - - - - - - - - - - - - - - - -
    -
    -
    - - -
    - -
    -

    src package

    -
    -

    Submodules

    -
    -
    -

    src.data module

    -
    -
    -

    src.date module

    -
    -
    -

    src.fit module

    -
    -
    -

    src.geometry module

    -
    -
    -

    Module contents

    -
    -
    - - -
    - -
    -
    - -
    -
    -
    - ©2019, Daniel Weschke. - - | - Powered by Sphinx 2.0.1 - & Alabaster 0.7.12 - - | - Page source -
    - - - - - - \ No newline at end of file diff --git a/docs/source/conf.py b/docs/source/conf.py index 6c43e5d..75aecc0 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -76,4 +76,6 @@ html_theme_options = { 'highlight_bg': '#444F65', 'xref_bg': 'transparent', 'xref_border': 'transparent', + 'seealso_bg': '#3c3c3c', + 'seealso_border': '#2C2C2C', } diff --git a/docs/source/modules.rst b/docs/source/modules.rst index 8b02d0a..6a908e5 100644 --- a/docs/source/modules.rst +++ b/docs/source/modules.rst @@ -8,3 +8,5 @@ src date fit geometry + solver + solver_model diff --git a/docs/source/solver.rst b/docs/source/solver.rst new file mode 100644 index 0000000..0d6d757 --- /dev/null +++ b/docs/source/solver.rst @@ -0,0 +1,7 @@ +solver module +============= + +.. automodule:: solver + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/solver_model.rst b/docs/source/solver_model.rst new file mode 100644 index 0000000..bc08685 --- /dev/null +++ b/docs/source/solver_model.rst @@ -0,0 +1,7 @@ +solver\_model module +==================== + +.. automodule:: solver_model + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/src.rst b/docs/source/src.rst deleted file mode 100644 index aa0f83a..0000000 --- a/docs/source/src.rst +++ /dev/null @@ -1,46 +0,0 @@ -src package -=========== - -Submodules ----------- - -src.data module ---------------- - -.. automodule:: src.data - :members: - :undoc-members: - :show-inheritance: - -src.date module ---------------- - -.. automodule:: src.date - :members: - :undoc-members: - :show-inheritance: - -src.fit module --------------- - -.. automodule:: src.fit - :members: - :undoc-members: - :show-inheritance: - -src.geometry module -------------------- - -.. automodule:: src.geometry - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: src - :members: - :undoc-members: - :show-inheritance: diff --git a/src/data.py b/src/data.py index 5f84c99..1c712d9 100644 --- a/src/data.py +++ b/src/data.py @@ -1,8 +1,29 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +"""Read and write data to or from file. + +.. module:: data + :platform: *nix, Windows + :synopsis: Handle data files. + +.. moduleauthor:: Daniel Weschke +""" 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()