daniel/pylib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Files
4fc4903dc204a295ba9abc863e2e14fd0a4d8d1f
pylib/docs/build/html/_modules/pylib/numerical/ode.html
Daniel Weschke 4fc4903dc2 change data, mathematics, function, geometry and geometry_plot_pylab. add data_step, data_step_std and helper. add an example and more documentation. 
inside the data module rename read to read_columns add add new read function to read the whole file as string. add print_list function to print one element per line. add unique_list and unique_list_hashable to reduce the list into a unique list with same order. add find_last, str_between, str_to_list functions.

inside the mathematics module for vector add normalized creating a new object (normalize will change the object), isclose and iscloseto, change ang to round internal number. for matrix improve slicing and add transposed creating a new object (transpose will change object).

inside the function module add b_spline_basis, b_spline_curve_with_knots and b_spline_knots functions. add sample_hal_open and sample_half_open_seq. add circle and ellipse.

inside the geometry module change CS init from using lists to Directions and add new constructor CS.init_xzy using lists. rename Wireframe to Polyline. add classes B_spline_curve_with_knots, Ellipse, ArcCircle, ArcEllipse, ArcBSplineCurveWithKnots. add function sample_half_open to create a list of Points.

inside the geometry_plot_pylab module change the help text.

add step_and data data_step_std module to read a step file to list and draw the content.

add helper module with context manager and decorator timeit to meassure the time for a section or function.

add example for b_spline function.
2020-01-08 21:59:53 +01:00

588 lines
64 KiB
HTML
Raw Blame History

<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta charset="utf-8" />
<title>pylib.numerical.ode &#8212; pylib 2019.12.21 documentation</title>
<link rel="stylesheet" href="../../../_static/alabaster.css" type="text/css" />
<link rel="stylesheet" href="../../../_static/pygments.css" type="text/css" />
<link rel="stylesheet" href="../../../_static/custom.css" type="text/css" />
<script type="text/javascript" id="documentation_options" data-url_root="../../../" src="../../../_static/documentation_options.js"></script>
<script type="text/javascript" src="../../../_static/jquery.js"></script>
<script type="text/javascript" src="../../../_static/underscore.js"></script>
<script type="text/javascript" src="../../../_static/doctools.js"></script>
<script type="text/javascript" src="../../../_static/language_data.js"></script>
<script async="async" type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
<link rel="index" title="Index" href="../../../genindex.html" />
<link rel="search" title="Search" href="../../../search.html" />
<link rel="stylesheet" href="../../../_static/custom.css" type="text/css" />
<meta name="viewport" content="width=device-width, initial-scale=0.9, maximum-scale=0.9" />
</head><body>
<div class="document">
<div class="documentwrapper">
<div class="bodywrapper">
<div class="body" role="main">
<h1>Source code for pylib.numerical.ode</h1><div class="highlight"><pre>
<span></span><span class="ch">#!/usr/bin/env python</span>
<span class="c1"># -*- coding: utf-8 -*-</span>
<span class="sa">r</span><span class="sd">&quot;&quot;&quot;Numerical solver of ordinary differential equations.</span>
<span class="sd">Solves the initial value problem for systems of first order</span>
<span class="sd">ordinary differential equations.</span>
<span class="sd">:Date: 2015-09-21</span>
<span class="sd">.. module:: ode</span>
<span class="sd"> :platform: *nix, Windows</span>
<span class="sd"> :synopsis: Numerical solver.</span>
<span class="sd">.. moduleauthor:: Daniel Weschke &lt;daniel.weschke@directbox.de&gt;</span>
<span class="sd">Approximate the solution :math:`x(t)` of the initial value problem</span>
<span class="sd">.. math ::</span>
<span class="sd"> \frac{\mathrm{d}x}{\mathrm{d}t} = \dot{x} &amp;= f(t,x) \\</span>
<span class="sd"> x(t_0) &amp;= x_0 \\</span>
<span class="sd"> t &amp;\in [t_0, t_n]</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="kn">from</span> <span class="nn">numpy</span> <span class="kn">import</span> <span class="n">array</span><span class="p">,</span> <span class="n">isnan</span><span class="p">,</span> <span class="nb">sum</span><span class="p">,</span> <span class="n">zeros</span><span class="p">,</span> <span class="n">dot</span>
<span class="kn">from</span> <span class="nn">numpy.linalg</span> <span class="kn">import</span> <span class="n">norm</span><span class="p">,</span> <span class="n">inv</span>
<div class="viewcode-block" id="e1"><a class="viewcode-back" href="../../../pylib.numerical.ode.html#pylib.numerical.ode.e1">[docs]</a><span class="k">def</span> <span class="nf">e1</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">x0</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
<span class="sa">r</span><span class="sd">&quot;&quot;&quot;Explicit first-order method /</span>
<span class="sd"> (standard, or forward) Euler method /</span>
<span class="sd"> Runge-Kutta 1st order method.</span>
<span class="sd"> de:</span>
<span class="sd"> Euler&#39;sche Polygonzugverfahren / explizite Euler-Verfahren /</span>
<span class="sd"> Euler-Cauchy-Verfahren / Euler-vorwärts-Verfahren</span>
<span class="sd"> :param f: the function to solve</span>
<span class="sd"> :type f: function</span>
<span class="sd"> :param x0: initial condition</span>
<span class="sd"> :type x0: list</span>
<span class="sd"> :param t: time</span>
<span class="sd"> :type t: list</span>
<span class="sd"> :param `*p`: parameters of the function (thickness, diameter,</span>
<span class="sd"> ...)</span>
<span class="sd"> :param verbose: print information (default = False)</span>
<span class="sd"> :type verbose: bool</span>
<span class="sd"> Approximate the solution :math:`x(t)` of the initial value problem</span>
<span class="sd"> .. math ::</span>
<span class="sd"> \dot{x} &amp;= f(t,x) \\</span>
<span class="sd"> x(t_0) &amp;= x_0 \\</span>
<span class="sd"> t &amp;\in [t_0, t_n]</span>
<span class="sd"> Choose a value h for the size of every step and set</span>
<span class="sd"> .. math ::</span>
<span class="sd"> t_{i+1} = t_0 + i h = t_i + h ~,\quad i=0,1,2,\ldots,n-1</span>
<span class="sd"> The derivative of the solution is approximated as the forward</span>
<span class="sd"> difference equation</span>
<span class="sd"> .. math ::</span>
<span class="sd"> \dot{x}_i = f(t_i, x_i) = \frac{x_{i+1} - x_i}{t_{i+1}-t_i}</span>
<span class="sd"> Therefore one step :math:`h` of the Euler method from</span>
<span class="sd"> :math:`t_i` to :math:`t_{i+1}` is</span>
<span class="sd"> .. math ::</span>
<span class="sd"> x_{i+1} &amp;= x_i + (t_{i+1}-t_i) f(t_i, x_i) \\</span>
<span class="sd"> x_{i+1} &amp;= x_i + h f(t_i, x_i) \\</span>
<span class="sd"> Example 1:</span>
<span class="sd"> .. math ::</span>
<span class="sd"> m\ddot{u} + d\dot{u} + ku = f(t) \\</span>
<span class="sd"> \ddot{u} = m^{-1}(f(t) - d\dot{u} - ku) \\</span>
<span class="sd"> with</span>
<span class="sd"> .. math ::</span>
<span class="sd"> x_1 &amp;= u &amp;\quad \dot{x}_1 = \dot{u} = x_2 \\</span>
<span class="sd"> x_2 &amp;= \dot{u} &amp;\quad \dot{x}_2 = \ddot{u} \\</span>
<span class="sd"> becomes</span>
<span class="sd"> .. math ::</span>
<span class="sd"> \dot{x}_1 &amp;= x_2 \\</span>
<span class="sd"> \dot{x}_2 &amp;= m^{-1}(f(t) - d x_2 - k x_1) \\</span>
<span class="sd"> or</span>
<span class="sd"> .. math ::</span>
<span class="sd"> \dot{x} &amp;= f(t,x) \\</span>
<span class="sd"> \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} &amp;=</span>
<span class="sd"> \begin{bmatrix} x_2 \\ m^{-1}(f(t) - d x_2 - k x_1)</span>
<span class="sd"> \end{bmatrix} \\</span>
<span class="sd"> &amp;=</span>
<span class="sd"> \begin{bmatrix} 0 \\ m^{-1} f(t) \end{bmatrix} +</span>
<span class="sd"> \begin{bmatrix} 0 &amp; 1 \\ -m^{-1} k &amp; -m^{-1} d \end{bmatrix}</span>
<span class="sd"> \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}</span>
<span class="sd"> Example 2:</span>
<span class="sd"> .. math ::</span>
<span class="sd"> m(u)\ddot{u} + d(u,\dot{u})\dot{u} + k(u)u = f(t) \\</span>
<span class="sd"> \ddot{u} = m^{-1}(u)(f(t) - d(u,\dot{u})\dot{u} - k(u)u) \\</span>
<span class="sd"> with</span>
<span class="sd"> .. math ::</span>
<span class="sd"> x_1 &amp;= u &amp;\quad \dot{x}_1 = \dot{u} = x_2 \\</span>
<span class="sd"> x_2 &amp;= \dot{u} &amp;\quad \dot{x}_2 = \ddot{u} \\</span>
<span class="sd"> becomes</span>
<span class="sd"> .. math ::</span>
<span class="sd"> \dot{x}_1 &amp;= x_2 \\</span>
<span class="sd"> \dot{x}_2 &amp;=</span>
<span class="sd"> m^{-1}(x_1)(f(t) - d(x_1,x_2) x_2 - k(x_1) x_1) \\</span>
<span class="sd"> or</span>
<span class="sd"> .. math ::</span>
<span class="sd"> \dot{x} &amp;= f(t,x) \\</span>
<span class="sd"> \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} &amp;=</span>
<span class="sd"> \begin{bmatrix}</span>
<span class="sd"> x_2 \\ m^{-1}(x_1)(f(t) - d(x_1,x_2) x_2 - k(x_1) x_1)</span>
<span class="sd"> \end{bmatrix} \\</span>
<span class="sd"> &amp;=</span>
<span class="sd"> \begin{bmatrix} 0 \\ m^{-1}(x_1) f(t) \end{bmatrix} +</span>
<span class="sd"> \begin{bmatrix}</span>
<span class="sd"> 0 &amp; 1 \\ -m^{-1}(x_1) k(x_1) &amp; -m^{-1} d(x_1,x_2)</span>
<span class="sd"> \end{bmatrix}</span>
<span class="sd"> \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}</span>
<span class="sd"> The Euler method is a first-order method, which means that the</span>
<span class="sd"> local error (error per step) is proportional to the square of</span>
<span class="sd"> the step size, and the global error (error at a given time) is</span>
<span class="sd"> proportional to the step size.</span>
<span class="sd"> &quot;&quot;&quot;</span>
<span class="kn">from</span> <span class="nn">pylib.data</span> <span class="kn">import</span> <span class="n">issequence</span>
<span class="k">if</span> <span class="n">issequence</span><span class="p">(</span><span class="n">x0</span><span class="p">):</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">x0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">x</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x0</span> <span class="c1"># Initial condition gives solution at first t</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span> <span class="c1"># Calculation loop</span>
<span class="n">Dt</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
<span class="n">dxdt</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:],</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="o">*</span><span class="n">p</span><span class="p">))</span>
<span class="c1"># Approximate solution at next value of x</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span> <span class="o">+</span> <span class="n">dxdt</span><span class="o">*</span><span class="n">Dt</span>
<span class="k">else</span><span class="p">:</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">)))</span>
<span class="n">x</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">x0</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span> <span class="c1"># Calculation loop</span>
<span class="n">Dt</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
<span class="n">dxdt</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="o">*</span><span class="n">p</span><span class="p">))</span>
<span class="c1"># Approximate solution at next value of x</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="n">dxdt</span><span class="o">*</span><span class="n">Dt</span>
<span class="k">if</span> <span class="n">verbose</span><span class="p">:</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Numerical integration of ODE using explicit &#39;</span> <span class="o">+</span>
<span class="s1">&#39;first-order method (Euler / Runge-Kutta) was successful.&#39;</span><span class="p">)</span>
<span class="k">return</span> <span class="n">x</span></div>
<div class="viewcode-block" id="e2"><a class="viewcode-back" href="../../../pylib.numerical.ode.html#pylib.numerical.ode.e2">[docs]</a><span class="k">def</span> <span class="nf">e2</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">x0</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
<span class="sa">r</span><span class="sd">&quot;&quot;&quot;Explicit second-order method / Runge-Kutta 2nd order method.</span>
<span class="sd"> :param f: the function to solve</span>
<span class="sd"> :type f: function</span>
<span class="sd"> :param x0: initial condition</span>
<span class="sd"> :type x0: list</span>
<span class="sd"> :param t: time</span>
<span class="sd"> :type t: list</span>
<span class="sd"> :param `*p`: parameters of the function (thickness, diameter,</span>
<span class="sd"> ...)</span>
<span class="sd"> :param verbose: print information (default = False)</span>
<span class="sd"> :type verbose: bool</span>
<span class="sd"> &quot;&quot;&quot;</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">x0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">x</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x0</span> <span class="c1"># Initial condition gives solution at first t</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span> <span class="c1"># Calculation loop</span>
<span class="n">Dt</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
<span class="n">k_1</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:],</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="o">*</span><span class="n">p</span><span class="p">))</span>
<span class="n">k_2</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">Dt</span><span class="o">*</span><span class="n">k_1</span><span class="p">,</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">Dt</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">))</span>
<span class="c1"># Approximate solution at next value of x</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span> <span class="o">+</span> <span class="n">k_2</span><span class="o">*</span><span class="n">Dt</span>
<span class="k">if</span> <span class="n">verbose</span><span class="p">:</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Numerical integration of ODE using explicit &#39;</span> <span class="o">+</span>
<span class="s1">&#39;2th-order method (Runge-Kutta) was successful.&#39;</span><span class="p">)</span>
<span class="k">return</span> <span class="n">x</span></div>
<div class="viewcode-block" id="e4"><a class="viewcode-back" href="../../../pylib.numerical.ode.html#pylib.numerical.ode.e4">[docs]</a><span class="k">def</span> <span class="nf">e4</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">x0</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
<span class="sa">r</span><span class="sd">&quot;&quot;&quot;Explicit fourth-order method / Runge-Kutta 4th order method.</span>
<span class="sd"> :param f: the function to solve</span>
<span class="sd"> :type f: function</span>
<span class="sd"> :param x0: initial condition</span>
<span class="sd"> :type x0: list</span>
<span class="sd"> :param t: time</span>
<span class="sd"> :type t: list</span>
<span class="sd"> :param `*p`: parameters of the function (thickness, diameter,</span>
<span class="sd"> ...)</span>
<span class="sd"> :param verbose: print information (default = False)</span>
<span class="sd"> :type verbose: bool</span>
<span class="sd"> Problem</span>
<span class="sd"> .. math ::</span>
<span class="sd"> \dot{y} &amp;= f(t, y) \\</span>
<span class="sd"> y(t_0) &amp;= y_0 \\</span>
<span class="sd"> t &amp;\in [t_0, t_n]</span>
<span class="sd"> Increment :math:`\delta t = t_{i+1}-t_i ~,~~ i=0,1,2,\ldots,n-1`</span>
<span class="sd"> .. math::</span>
<span class="sd"> y_{n+1} &amp;= y_{i} + \tfrac{1}{6}(</span>
<span class="sd"> \delta y_{i,1} + 2\delta y_{i,2} + 2\delta y_{i,3} +</span>
<span class="sd"> \delta y_{i,4}) \\</span>
<span class="sd"> &amp; \qquad \text{with} \\</span>
<span class="sd"> \delta y_{i,1} &amp;= \delta t \cdot y&#39;(t_{i}, ~ y_{i}) \\</span>
<span class="sd"> \delta y_{i,2} &amp;= \delta t \cdot</span>
<span class="sd"> y&#39;(t_{i}+\tfrac{1}{2}\delta t, ~</span>
<span class="sd"> y_{i}+\tfrac{1}{2}\delta y_{i,1}) \\</span>
<span class="sd"> \delta y_{i,3} &amp;= \delta t \cdot</span>
<span class="sd"> y&#39;(t_{i}+\tfrac{1}{2}\delta t, ~</span>
<span class="sd"> y_{i}+\tfrac{1}{2}\delta y_{i,2}) \\</span>
<span class="sd"> \delta y_{i,4} &amp;= \delta t \cdot</span>
<span class="sd"> y&#39;(t_{i}+\delta t, ~ y_{i}+\delta y_{i,3})</span>
<span class="sd"> &quot;&quot;&quot;</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">x0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">x</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x0</span> <span class="c1"># Initial condition</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span> <span class="c1"># Calculation loop</span>
<span class="n">Dt</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
<span class="n">k_1</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:],</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="o">*</span><span class="n">p</span><span class="p">))</span>
<span class="n">k_2</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">Dt</span><span class="o">*</span><span class="n">k_1</span><span class="p">,</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">Dt</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">))</span>
<span class="n">k_3</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">Dt</span><span class="o">*</span><span class="n">k_2</span><span class="p">,</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">+</span><span class="mf">0.5</span><span class="o">*</span><span class="n">Dt</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">))</span>
<span class="n">k_4</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">f</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">+</span><span class="n">k_3</span><span class="o">*</span><span class="n">Dt</span><span class="p">,</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">+</span><span class="n">Dt</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">))</span>
<span class="c1"># Approximate solution at next value of x</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span> <span class="o">+</span> <span class="mf">1.</span><span class="o">/</span><span class="mi">6</span><span class="o">*</span><span class="p">(</span><span class="n">k_1</span><span class="o">+</span><span class="mi">2</span><span class="o">*</span><span class="n">k_2</span><span class="o">+</span><span class="mi">2</span><span class="o">*</span><span class="n">k_3</span><span class="o">+</span><span class="n">k_4</span><span class="p">)</span><span class="o">*</span><span class="n">Dt</span>
<span class="k">if</span> <span class="n">verbose</span><span class="p">:</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Numerical integration of ODE using explicit &#39;</span> <span class="o">+</span>
<span class="s1">&#39;4th-order method (Runge-Kutta) was successful.&#39;</span><span class="p">)</span>
<span class="k">return</span> <span class="n">x</span></div>
<div class="viewcode-block" id="fpi"><a class="viewcode-back" href="../../../pylib.numerical.ode.html#pylib.numerical.ode.fpi">[docs]</a><span class="k">def</span> <span class="nf">fpi</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">xi</span><span class="p">,</span> <span class="n">ti</span><span class="p">,</span> <span class="n">ti1</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">,</span> <span class="n">max_iterations</span><span class="o">=</span><span class="mi">1000</span><span class="p">,</span> <span class="n">tol</span><span class="o">=</span><span class="mf">1e-9</span><span class="p">,</span>
<span class="n">verbose</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
<span class="sa">r</span><span class="sd">&quot;&quot;&quot;Fixed-point iteration.</span>
<span class="sd"> :param f: the function to iterate :math:`f = \dot{x}(x,t)`</span>
<span class="sd"> :type f: function</span>
<span class="sd"> :param xi: initial condition :math:`x_i`</span>
<span class="sd"> :type xi: list</span>
<span class="sd"> :param ti: time :math:`t_i`</span>
<span class="sd"> :type ti: float</span>
<span class="sd"> :param ti1: time :math:`t_{i+1}`</span>
<span class="sd"> :type ti1: float</span>
<span class="sd"> :param `*p`: parameters of the function (thickness, diameter,</span>
<span class="sd"> ...)</span>
<span class="sd"> :param max_iterations: maximum number of iterations</span>
<span class="sd"> :type max_iterations: int</span>
<span class="sd"> :param tol: tolerance against residuum :math:`\varepsilon`</span>
<span class="sd"> (default = 1e-9)</span>
<span class="sd"> :type tol: float</span>
<span class="sd"> :param verbose: print information (default = False)</span>
<span class="sd"> :type verbose: bool</span>
<span class="sd"> :returns: :math:`x_{i}`</span>
<span class="sd"> .. math ::</span>
<span class="sd"> x_{i,j=0} = x_{i}</span>
<span class="sd"> .. math ::</span>
<span class="sd"> x_{i,j+1} = x_i + \dot{x}(x_{i,j}, t_{i+1})\cdot(t_{i+1}-t_i)</span>
<span class="sd"> .. math ::</span>
<span class="sd"> \text{residuum} = \frac{\lVert x_{i,j+1}-x_{i,j}\rVert}</span>
<span class="sd"> {\lVert x_{i,j+1} \rVert} &lt; \varepsilon</span>
<span class="sd"> .. math ::</span>
<span class="sd"> x_{i} = x_{i,j=\text{end}}</span>
<span class="sd"> &quot;&quot;&quot;</span>
<span class="n">xij</span> <span class="o">=</span> <span class="n">xi</span>
<span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">max_iterations</span><span class="p">):</span>
<span class="n">dxdt</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">f</span><span class="p">(</span><span class="n">xij</span><span class="p">,</span> <span class="n">ti1</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">))</span>
<span class="c1"># Approximate solution at next value of x</span>
<span class="n">xij1</span> <span class="o">=</span> <span class="n">xi</span> <span class="o">+</span> <span class="n">dxdt</span> <span class="o">*</span> <span class="p">(</span><span class="n">ti1</span><span class="o">-</span><span class="n">ti</span><span class="p">)</span>
<span class="n">residuum</span> <span class="o">=</span> <span class="n">norm</span><span class="p">(</span><span class="n">xij1</span><span class="o">-</span><span class="n">xij</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">xij1</span><span class="p">)</span>
<span class="n">xij</span> <span class="o">=</span> <span class="n">xij1</span>
<span class="k">if</span> <span class="n">residuum</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">:</span>
<span class="k">break</span>
<span class="n">iterations</span> <span class="o">=</span> <span class="n">j</span><span class="o">+</span><span class="mi">1</span> <span class="c1"># number beginning with 1 therefore + 1</span>
<span class="k">return</span> <span class="n">xij</span><span class="p">,</span> <span class="n">iterations</span></div>
<div class="viewcode-block" id="i1"><a class="viewcode-back" href="../../../pylib.numerical.ode.html#pylib.numerical.ode.i1">[docs]</a><span class="k">def</span> <span class="nf">i1</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">x0</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">,</span> <span class="n">max_iterations</span><span class="o">=</span><span class="mi">1000</span><span class="p">,</span> <span class="n">tol</span><span class="o">=</span><span class="mf">1e-9</span><span class="p">,</span>
<span class="n">verbose</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
<span class="sa">r</span><span class="sd">&quot;&quot;&quot;Implicite first-order method / backward Euler method.</span>
<span class="sd"> :param f: the function to solve</span>
<span class="sd"> :type f: function</span>
<span class="sd"> :param x0: initial condition</span>
<span class="sd"> :type x0: list</span>
<span class="sd"> :param t: time</span>
<span class="sd"> :type t: list</span>
<span class="sd"> :param `*p`: parameters of the function (thickness, diameter,</span>
<span class="sd"> ...)</span>
<span class="sd"> :param max_iterations: maximum number of iterations</span>
<span class="sd"> :type max_iterations: int</span>
<span class="sd"> :param tol: tolerance against residuum (default = 1e-9)</span>
<span class="sd"> :type tol: float</span>
<span class="sd"> :param verbose: print information (default = False)</span>
<span class="sd"> :type verbose: bool</span>
<span class="sd"> The backward Euler method has order one and is A-stable.</span>
<span class="sd"> &quot;&quot;&quot;</span>
<span class="n">iterations</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="mi">1</span><span class="p">))</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">x0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">x</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x0</span> <span class="c1"># Initial condition gives solution at first t</span>
<span class="c1"># x(i+1) = x(i) + f(x(i+1), t(i+1)), exact value of</span>
<span class="c1"># f(x(i+1), t(i+1)) is not available therefore using</span>
<span class="c1"># Newton-Raphson method</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
<span class="n">Dt</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
<span class="n">xi</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span>
<span class="n">xi</span><span class="p">,</span> <span class="n">iteration</span> <span class="o">=</span> <span class="n">fpi</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">xi</span><span class="p">,</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">],</span> <span class="o">*</span><span class="n">p</span><span class="p">,</span>
<span class="n">max_iterations</span><span class="o">=</span><span class="n">max_iterations</span><span class="p">,</span> <span class="n">tol</span><span class="o">=</span><span class="n">tol</span><span class="p">,</span>
<span class="n">verbose</span><span class="o">=</span><span class="n">verbose</span><span class="p">)</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">xi</span>
<span class="n">iterations</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">iteration</span>
<span class="k">if</span> <span class="n">verbose</span><span class="p">:</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Numerical integration of ODE using implicite &#39;</span> <span class="o">+</span>
<span class="s1">&#39;first-order method (Euler) was successful.&#39;</span><span class="p">)</span>
<span class="k">return</span> <span class="n">x</span><span class="p">,</span> <span class="n">iterations</span></div>
<div class="viewcode-block" id="newmark_newtonraphson"><a class="viewcode-back" href="../../../pylib.numerical.ode.html#pylib.numerical.ode.newmark_newtonraphson">[docs]</a><span class="k">def</span> <span class="nf">newmark_newtonraphson</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">x0</span><span class="p">,</span> <span class="n">xp0</span><span class="p">,</span> <span class="n">xpp0</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">,</span> <span class="n">gamma</span><span class="o">=.</span><span class="mi">5</span><span class="p">,</span>
<span class="n">beta</span><span class="o">=.</span><span class="mi">25</span><span class="p">,</span> <span class="n">max_iterations</span><span class="o">=</span><span class="mi">1000</span><span class="p">,</span> <span class="n">tol</span><span class="o">=</span><span class="mf">1e-9</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
<span class="sa">r</span><span class="sd">&quot;&quot;&quot;Newmark method.</span>
<span class="sd"> :param f: the function to solve</span>
<span class="sd"> :type f: function</span>
<span class="sd"> :param x0: initial condition</span>
<span class="sd"> :type x0: list</span>
<span class="sd"> :param xp0: initial condition</span>
<span class="sd"> :type xp0: list</span>
<span class="sd"> :param xpp0: initial condition</span>
<span class="sd"> :type xpp0: list</span>
<span class="sd"> :param t: time</span>
<span class="sd"> :type t: list</span>
<span class="sd"> :param `*p`: parameters of the function (thickness, diameter,</span>
<span class="sd"> ...)</span>
<span class="sd"> :param gamma: newmark parameter for velocity (default = 0.5)</span>
<span class="sd"> :type gamma: float</span>
<span class="sd"> :param beta: newmark parameter for displacement (default = 0.25)</span>
<span class="sd"> :type beta: float</span>
<span class="sd"> :param max_iterations: maximum number of iterations</span>
<span class="sd"> :type max_iterations: int</span>
<span class="sd"> :param tol: tolerance against residuum (default = 1e-9)</span>
<span class="sd"> :type tol: float</span>
<span class="sd"> :param verbose: print information (default = False)</span>
<span class="sd"> :type verbose: bool</span>
<span class="sd"> &quot;&quot;&quot;</span>
<span class="n">iterations</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="mi">1</span><span class="p">))</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">x0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">xp</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">xp0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">xpp</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">xpp0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">x</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x0</span> <span class="c1"># Initial condition gives solution at first t</span>
<span class="n">xp</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">xp0</span> <span class="c1"># Initial condition gives solution at first t</span>
<span class="n">xpp</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">xpp0</span> <span class="c1"># Initial condition gives solution at first t</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
<span class="n">Dt</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
<span class="n">xi</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">xpi</span> <span class="o">=</span> <span class="n">xp</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">xppi</span> <span class="o">=</span> <span class="n">xpp</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">x1</span> <span class="o">=</span> <span class="n">xi</span>
<span class="n">xp1</span> <span class="o">=</span> <span class="n">xpi</span>
<span class="n">xpp1</span> <span class="o">=</span> <span class="n">xppi</span>
<span class="n">j</span> <span class="o">=</span> <span class="mi">0</span>
<span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">max_iterations</span><span class="p">):</span> <span class="c1"># Fixed-point iteration</span>
<span class="c1">#dxdt = array(f(t[i+1], x1, p))</span>
<span class="c1"># Approximate solution at next value of x</span>
<span class="c1">#x11 = x[i,:] + dxdt*Dt</span>
<span class="n">N</span><span class="p">,</span> <span class="n">dN</span><span class="p">,</span> <span class="n">dNp</span><span class="p">,</span> <span class="n">dNpp</span> <span class="o">=</span> <span class="n">f</span><span class="p">(</span><span class="n">x1</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,)</span><span class="o">.</span><span class="n">tolist</span><span class="p">(),</span>
<span class="n">xp1</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,)</span><span class="o">.</span><span class="n">tolist</span><span class="p">(),</span> <span class="n">xpp1</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,)</span><span class="o">.</span><span class="n">tolist</span><span class="p">(),</span>
<span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="o">*</span><span class="n">p</span><span class="p">)</span>
<span class="k">if</span> <span class="n">isnan</span><span class="p">(</span><span class="nb">sum</span><span class="p">(</span><span class="n">dN</span><span class="p">))</span> <span class="ow">or</span> <span class="n">isnan</span><span class="p">(</span><span class="nb">sum</span><span class="p">(</span><span class="n">dNp</span><span class="p">))</span> <span class="ow">or</span> <span class="n">isnan</span><span class="p">(</span><span class="nb">sum</span><span class="p">(</span><span class="n">dNpp</span><span class="p">)):</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;divergiert&#39;</span><span class="p">)</span>
<span class="k">break</span>
<span class="n">xpp11</span> <span class="o">=</span> <span class="n">xpp1</span> <span class="o">-</span> <span class="n">dot</span><span class="p">(</span><span class="n">inv</span><span class="p">(</span><span class="n">dNpp</span><span class="p">),</span> <span class="p">(</span><span class="n">N</span> <span class="o">+</span> <span class="n">dot</span><span class="p">(</span><span class="n">dN</span><span class="p">,</span> <span class="p">(</span><span class="n">x1</span><span class="o">-</span><span class="n">xi</span><span class="p">))</span> <span class="o">+</span> \
<span class="n">dot</span><span class="p">(</span><span class="n">dNp</span><span class="p">,</span> <span class="p">(</span><span class="n">xp1</span><span class="o">-</span><span class="n">xpi</span><span class="p">))))</span>
<span class="n">xp1</span> <span class="o">=</span> <span class="n">xpi</span> <span class="o">+</span> <span class="n">Dt</span><span class="o">*</span><span class="p">(</span> <span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">gamma</span><span class="p">)</span><span class="o">*</span><span class="n">xppi</span> <span class="o">+</span> <span class="n">gamma</span><span class="o">*</span><span class="n">xpp11</span> <span class="p">)</span>
<span class="n">x1</span> <span class="o">=</span> <span class="n">xi</span> <span class="o">+</span> <span class="n">Dt</span><span class="o">*</span><span class="n">xpi</span> <span class="o">+</span> <span class="n">Dt</span><span class="o">**</span><span class="mi">2</span><span class="o">*</span><span class="p">(</span> <span class="p">(</span><span class="o">.</span><span class="mi">5</span><span class="o">-</span><span class="n">beta</span><span class="p">)</span><span class="o">*</span><span class="n">xppi</span> <span class="o">+</span> <span class="n">beta</span><span class="o">*</span><span class="n">xpp11</span> <span class="p">)</span>
<span class="n">residuum</span> <span class="o">=</span> <span class="n">norm</span><span class="p">(</span><span class="n">xpp11</span><span class="o">-</span><span class="n">xpp1</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">xpp11</span><span class="p">)</span>
<span class="n">xpp1</span> <span class="o">=</span> <span class="n">xpp11</span>
<span class="k">if</span> <span class="n">residuum</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">:</span>
<span class="k">break</span>
<span class="n">iterations</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">j</span><span class="o">+</span><span class="mi">1</span>
<span class="n">xpp</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">xpp1</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,)</span><span class="o">.</span><span class="n">tolist</span><span class="p">()</span>
<span class="n">xp</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">xp1</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,)</span><span class="o">.</span><span class="n">tolist</span><span class="p">()</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x1</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,)</span><span class="o">.</span><span class="n">tolist</span><span class="p">()</span>
<span class="k">if</span> <span class="n">verbose</span><span class="p">:</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Numerical integration of ODE using explicite &#39;</span> <span class="o">+</span>
<span class="s1">&#39;newmark method was successful.&#39;</span><span class="p">)</span>
<span class="k">return</span> <span class="n">x</span><span class="p">,</span> <span class="n">xp</span><span class="p">,</span> <span class="n">xpp</span><span class="p">,</span> <span class="n">iterations</span></div>
<span class="c1"># x = concatenate((x, xp, xpp), axis=1)</span>
<div class="viewcode-block" id="newmark_newtonraphson_mdk"><a class="viewcode-back" href="../../../pylib.numerical.ode.html#pylib.numerical.ode.newmark_newtonraphson_mdk">[docs]</a><span class="k">def</span> <span class="nf">newmark_newtonraphson_mdk</span><span class="p">(</span><span class="n">fmdk</span><span class="p">,</span> <span class="n">x0</span><span class="p">,</span> <span class="n">xp0</span><span class="p">,</span> <span class="n">xpp0</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="o">*</span><span class="n">p</span><span class="p">,</span> <span class="n">gamma</span><span class="o">=.</span><span class="mi">5</span><span class="p">,</span>
<span class="n">beta</span><span class="o">=.</span><span class="mi">25</span><span class="p">,</span> <span class="n">max_iterations</span><span class="o">=</span><span class="mi">1000</span><span class="p">,</span> <span class="n">tol</span><span class="o">=</span><span class="mf">1e-9</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
<span class="sa">r</span><span class="sd">&quot;&quot;&quot;Newmark method.</span>
<span class="sd"> Using m mass, d damping and k stiffness formulation.</span>
<span class="sd"> :param f: the function to solve</span>
<span class="sd"> :type f: function</span>
<span class="sd"> :param x0: initial condition</span>
<span class="sd"> :type x0: list</span>
<span class="sd"> :param xp0: initial condition</span>
<span class="sd"> :type xp0: list</span>
<span class="sd"> :param xpp0: initial condition</span>
<span class="sd"> :type xpp0: list</span>
<span class="sd"> :param t: time</span>
<span class="sd"> :type t: list</span>
<span class="sd"> :param `*p`: parameters of the function (thickness, diameter,</span>
<span class="sd"> ...)</span>
<span class="sd"> :param gamma: newmark parameter for velocity (default = 0.5)</span>
<span class="sd"> :type gamma: float</span>
<span class="sd"> :param beta: newmark parameter for displacement (default = 0.25)</span>
<span class="sd"> :type beta: float</span>
<span class="sd"> :param max_iterations: maximum number of iterations</span>
<span class="sd"> :type max_iterations: int</span>
<span class="sd"> :param tol: tolerance against residuum (default = 1e-9)</span>
<span class="sd"> :type tol: float</span>
<span class="sd"> :param verbose: print information (default = False)</span>
<span class="sd"> :type verbose: bool</span>
<span class="sd"> &quot;&quot;&quot;</span>
<span class="n">iterations</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="mi">1</span><span class="p">))</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">x0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">xp</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">xp0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">xpp</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">),</span> <span class="nb">len</span><span class="p">(</span><span class="n">xpp0</span><span class="p">)))</span> <span class="c1"># Preallocate array</span>
<span class="n">x</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x0</span> <span class="c1"># Initial condition gives solution at first t</span>
<span class="n">xp</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">xp0</span> <span class="c1"># Initial condition gives solution at first t</span>
<span class="n">xpp</span><span class="p">[</span><span class="mi">0</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">xpp0</span> <span class="c1"># Initial condition gives solution at first t</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">t</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
<span class="n">Dt</span> <span class="o">=</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
<span class="n">rm</span><span class="p">,</span> <span class="n">rmx</span><span class="p">,</span> <span class="n">rmxpp</span><span class="p">,</span> <span class="n">rd</span><span class="p">,</span> <span class="n">rdx</span><span class="p">,</span> <span class="n">rdxp</span><span class="p">,</span> <span class="n">rk</span><span class="p">,</span> <span class="n">rkx</span><span class="p">,</span> <span class="n">f</span> <span class="o">=</span> <span class="n">fmdk</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:],</span>
<span class="n">xp</span><span class="p">[</span><span class="n">i</span><span class="p">,:],</span> <span class="n">xpp</span><span class="p">[</span><span class="n">i</span><span class="p">,:],</span> <span class="n">t</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="o">*</span><span class="n">p</span><span class="p">)</span>
<span class="n">xi</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">xpi</span> <span class="o">=</span> <span class="n">xp</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">xppi</span> <span class="o">=</span> <span class="n">xpp</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)</span>
<span class="n">x1</span> <span class="o">=</span> <span class="n">xi</span>
<span class="n">xp1</span> <span class="o">=</span> <span class="n">xpi</span>
<span class="n">xpp1</span> <span class="o">=</span> <span class="n">xppi</span>
<span class="n">j</span> <span class="o">=</span> <span class="mi">0</span>
<span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">max_iterations</span><span class="p">):</span> <span class="c1"># Fixed-point iteration</span>
<span class="c1">#dxdt = array(f(t[i+1], x1, p))</span>
<span class="c1"># Approximate solution at next value of x</span>
<span class="c1">#x11 = x[i,:] + dxdt*Dt</span>
<span class="n">r</span> <span class="o">=</span> <span class="p">(</span><span class="n">rmx</span><span class="o">+</span><span class="n">rdx</span><span class="o">+</span><span class="n">rkx</span><span class="p">)</span><span class="o">*</span><span class="n">Dt</span><span class="o">**</span><span class="mf">2.</span><span class="o">/</span><span class="mi">4</span> <span class="o">+</span> <span class="n">rdxp</span><span class="o">*</span><span class="n">Dt</span><span class="o">/</span><span class="mi">2</span> <span class="o">+</span> <span class="n">rmxpp</span>
<span class="n">rp</span> <span class="o">=</span> <span class="n">f</span> <span class="o">-</span> <span class="p">(</span><span class="n">rm</span> <span class="o">+</span> <span class="n">dot</span><span class="p">(</span><span class="n">rmx</span><span class="p">,</span> <span class="p">(</span><span class="n">Dt</span><span class="o">*</span><span class="n">xpi</span><span class="o">+</span><span class="n">Dt</span><span class="o">**</span><span class="mf">2.</span><span class="o">/</span><span class="mi">4</span><span class="o">*</span><span class="n">xppi</span><span class="p">))</span> <span class="o">-</span> \
<span class="n">dot</span><span class="p">(</span><span class="n">rmxpp</span><span class="p">,</span> <span class="n">xppi</span><span class="p">)</span> <span class="o">+</span> \
<span class="n">rd</span> <span class="o">+</span> <span class="n">dot</span><span class="p">(</span><span class="n">rdx</span><span class="p">,</span> <span class="p">(</span><span class="n">Dt</span><span class="o">*</span><span class="n">xpi</span><span class="o">+</span><span class="n">Dt</span><span class="o">**</span><span class="mf">2.</span><span class="o">/</span><span class="mi">4</span><span class="o">*</span><span class="n">xppi</span><span class="p">))</span> <span class="o">+</span> \
<span class="n">dot</span><span class="p">(</span><span class="n">rdxp</span><span class="p">,</span> <span class="n">Dt</span><span class="o">/</span><span class="mi">2</span><span class="o">*</span><span class="n">xppi</span><span class="p">)</span> <span class="o">+</span> \
<span class="n">rk</span> <span class="o">+</span> <span class="n">dot</span><span class="p">(</span><span class="n">rkx</span><span class="p">,</span> <span class="p">(</span><span class="n">Dt</span><span class="o">*</span><span class="n">xpi</span><span class="o">+</span><span class="n">Dt</span><span class="o">**</span><span class="mf">2.</span><span class="o">/</span><span class="mi">4</span><span class="o">*</span><span class="n">xppi</span><span class="p">))</span> <span class="p">)</span>
<span class="n">xpp11</span> <span class="o">=</span> <span class="n">dot</span><span class="p">(</span><span class="n">inv</span><span class="p">(</span><span class="n">r</span><span class="p">),</span> <span class="n">rp</span><span class="p">)</span>
<span class="n">xp1</span> <span class="o">=</span> <span class="n">xpi</span> <span class="o">+</span> <span class="n">Dt</span><span class="o">*</span><span class="p">(</span> <span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">gamma</span><span class="p">)</span><span class="o">*</span><span class="n">xppi</span> <span class="o">+</span> <span class="n">gamma</span><span class="o">*</span><span class="n">xpp11</span> <span class="p">)</span>
<span class="n">x1</span> <span class="o">=</span> <span class="n">xi</span> <span class="o">+</span> <span class="n">Dt</span><span class="o">*</span><span class="n">xpi</span> <span class="o">+</span> <span class="n">Dt</span><span class="o">**</span><span class="mi">2</span><span class="o">*</span><span class="p">(</span> <span class="p">(</span><span class="o">.</span><span class="mi">5</span><span class="o">-</span><span class="n">beta</span><span class="p">)</span><span class="o">*</span><span class="n">xppi</span> <span class="o">+</span> <span class="n">beta</span><span class="o">*</span><span class="n">xpp11</span> <span class="p">)</span>
<span class="n">residuum</span> <span class="o">=</span> <span class="n">norm</span><span class="p">(</span><span class="n">xpp11</span><span class="o">-</span><span class="n">xpp1</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">xpp11</span><span class="p">)</span>
<span class="n">xpp1</span> <span class="o">=</span> <span class="n">xpp11</span>
<span class="k">if</span> <span class="n">residuum</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">:</span>
<span class="k">break</span>
<span class="n">iterations</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">j</span><span class="o">+</span><span class="mi">1</span>
<span class="n">xpp</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">xpp1</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,)</span><span class="o">.</span><span class="n">tolist</span><span class="p">()</span>
<span class="n">xp</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">xp1</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,)</span><span class="o">.</span><span class="n">tolist</span><span class="p">()</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">x1</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,)</span><span class="o">.</span><span class="n">tolist</span><span class="p">()</span>
<span class="k">if</span> <span class="n">verbose</span><span class="p">:</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Numerical integration of ODE using explicite &#39;</span> <span class="o">+</span>
<span class="s1">&#39;newmark method was successful.&#39;</span><span class="p">)</span>
<span class="k">return</span> <span class="n">x</span><span class="p">,</span> <span class="n">xp</span><span class="p">,</span> <span class="n">xpp</span><span class="p">,</span> <span class="n">iterations</span></div>
<span class="c1"># x = concatenate((x, xp, xpp), axis=1)</span>
</pre></div>
</div>
</div>
</div>
<div class="sphinxsidebar" role="navigation" aria-label="main navigation">
<div class="sphinxsidebarwrapper">
<h1 class="logo"><a href="../../../index.html">pylib</a></h1>
<h3>Navigation</h3>
<p class="caption"><span class="caption-text">Contents:</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="../../../modules.html">pylib</a></li>
</ul>
<div class="relations">
<h3>Related Topics</h3>
<ul>
<li><a href="../../../index.html">Documentation overview</a><ul>
<li><a href="../../index.html">Module code</a><ul>
</ul></li>
</ul></li>
</ul>
</div>
<div id="searchbox" style="display: none" role="search">
<h3 id="searchlabel">Quick search</h3>
<div class="searchformwrapper">
<form class="search" action="../../../search.html" method="get">
<input type="text" name="q" aria-labelledby="searchlabel" />
<input type="submit" value="Go" />
</form>
</div>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
</div>
</div>
<div class="clearer"></div>
</div>
<div class="footer">
&copy;2019, Daniel Weschke.
|
Powered by <a href="http://sphinx-doc.org/">Sphinx 2.2.1</a>
&amp; <a href="https://github.com/bitprophet/alabaster">Alabaster 0.7.12</a>
</div>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink