add time module and move fixed-point iteration to own function
This commit is contained in:
221
docs/build/html/_modules/numerical/ode.html
vendored
221
docs/build/html/_modules/numerical/ode.html
vendored
@@ -37,8 +37,8 @@
|
||||
<span class="c1"># -*- coding: utf-8 -*-</span>
|
||||
<span class="sd">"""Numerical solver of ordinary differential equations.</span>
|
||||
|
||||
<span class="sd">Solves the initial value problem for systems of first order ordinary differential</span>
|
||||
<span class="sd">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>
|
||||
|
||||
@@ -68,7 +68,8 @@
|
||||
<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"> :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>
|
||||
|
||||
@@ -83,14 +84,14 @@
|
||||
<span class="sd"> .. math ::</span>
|
||||
<span class="sd"> t_i = t_0 + i h ~,\quad i=1,2,\ldots,n</span>
|
||||
|
||||
<span class="sd"> The derivative of the solution is approximated as the forward difference</span>
|
||||
<span class="sd"> equation</span>
|
||||
<span class="sd"> The derivative of the solution is approximated as the forward</span>
|
||||
<span class="sd"> difference equation</span>
|
||||
<span class="sd"> </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 :math:`t_i` to</span>
|
||||
<span class="sd"> :math:`t_{i+1}` is</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} &= x_i + (t_{i+1}-t_i) f(t_i, x_i) \\</span>
|
||||
@@ -119,7 +120,8 @@
|
||||
<span class="sd"> .. math ::</span>
|
||||
<span class="sd"> \dot{x} &= f(t,x) \\</span>
|
||||
<span class="sd"> \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} &=</span>
|
||||
<span class="sd"> \begin{bmatrix} x_2 \\ m^{-1}(f(t) - d x_2 - k x_1) \end{bmatrix} \\</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"> &=</span>
|
||||
<span class="sd"> \begin{bmatrix} 0 \\ m^{-1} f(t) \end{bmatrix} +</span>
|
||||
<span class="sd"> \begin{bmatrix} 0 & 1 \\ -m^{-1} k & -m^{-1} d \end{bmatrix}</span>
|
||||
@@ -141,32 +143,39 @@
|
||||
|
||||
<span class="sd"> .. math ::</span>
|
||||
<span class="sd"> \dot{x}_1 &= x_2 \\</span>
|
||||
<span class="sd"> \dot{x}_2 &= m^{-1}(x_1)(f(t) - d(x_1,x_2) x_2 - k(x_1) x_1) \\</span>
|
||||
<span class="sd"> \dot{x}_2 &=</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} &= f(t,x) \\</span>
|
||||
<span class="sd"> \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} &=</span>
|
||||
<span class="sd"> \begin{bmatrix} x_2 \\ m^{-1}(x_1)(f(t) - d(x_1,x_2) x_2 - k(x_1) x_1) \end{bmatrix} \\</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"> &=</span>
|
||||
<span class="sd"> \begin{bmatrix} 0 \\ m^{-1}(x_1) f(t) \end{bmatrix} +</span>
|
||||
<span class="sd"> \begin{bmatrix} 0 & 1 \\ -m^{-1}(x_1) k(x_1) & -m^{-1} d(x_1,x_2) \end{bmatrix}</span>
|
||||
<span class="sd"> \begin{bmatrix}</span>
|
||||
<span class="sd"> 0 & 1 \\ -m^{-1}(x_1) k(x_1) & -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,</span>
|
||||
<span class="sd"> which means that the local error (error per step) is proportional to the</span>
|
||||
<span class="sd"> square of the step size, and the global error (error at a given time) is</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"> """</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">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="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="c1"># Approximate solution at next value of x</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">'Numerical integration of ODE using explicit first-order method (Euler / Runge-Kutta) was successful.'</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Numerical integration of ODE using explicit '</span> <span class="o">+</span>
|
||||
<span class="s1">'first-order method (Euler / Runge-Kutta) was successful.'</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="../../numerical.html#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>
|
||||
@@ -178,19 +187,22 @@
|
||||
<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"> :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"> """</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">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="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="c1"># Approximate solution at next value of x</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">'Numerical integration of ODE using explicit 2th-order method (Runge-Kutta) was successful.'</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Numerical integration of ODE using explicit '</span> <span class="o">+</span>
|
||||
<span class="s1">'2th-order method (Runge-Kutta) was successful.'</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="../../numerical.html#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>
|
||||
@@ -202,7 +214,8 @@
|
||||
<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"> :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"> """</span>
|
||||
@@ -214,70 +227,64 @@
|
||||
<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="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="c1"># Approximate solution at next value of x</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">'Numerical integration of ODE using explicit 4th-order method (Runge-Kutta) was successful.'</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Numerical integration of ODE using explicit '</span> <span class="o">+</span>
|
||||
<span class="s1">'4th-order method (Runge-Kutta) was successful.'</span><span class="p">)</span>
|
||||
<span class="k">return</span> <span class="n">x</span></div>
|
||||
|
||||
<div class="viewcode-block" id="dxdt_Dt"><a class="viewcode-back" href="../../numerical.html#numerical.ode.dxdt_Dt">[docs]</a><span class="k">def</span> <span class="nf">dxdt_Dt</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">t</span><span class="p">,</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="sa">r</span><span class="sd">"""</span>
|
||||
<span class="sd"> :param f: :math:`f = \dot{x}`</span>
|
||||
<span class="sd"> :type f: function</span>
|
||||
<span class="sd"> :param Dt: :math:`\Delta{t}`</span>
|
||||
<div class="viewcode-block" id="fpi"><a class="viewcode-back" href="../../numerical.html#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">"""Fixed-point iteration.</span>
|
||||
<span class="sd"> </span>
|
||||
<span class="sd"> :returns: :math:`\Delta x = \dot{x} \Delta t`</span>
|
||||
<span class="sd"> """</span>
|
||||
<span class="k">return</span> <span class="n">array</span><span class="p">(</span><span class="n">dxdt</span><span class="p">(</span><span class="n">x</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="o">*</span> <span class="n">Dt</span></div>
|
||||
|
||||
<div class="viewcode-block" id="fixed_point_iteration"><a class="viewcode-back" href="../../numerical.html#numerical.ode.fixed_point_iteration">[docs]</a><span class="k">def</span> <span class="nf">fixed_point_iteration</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">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">"""</span>
|
||||
<span class="sd"> :param f: the function to iterate :math:`f = \Delta{x}(t)`</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 t: time :math:`t`</span>
|
||||
<span class="sd"> :type t: float</span>
|
||||
<span class="sd"> :param `*p`: parameters of the function (thickness, diameter, ...)</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 (default = 1e-9)</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"> </span>
|
||||
<span class="sd"> :returns: :math:`x_{i+1}`</span>
|
||||
<span class="sd"> :returns: :math:`x_{i}`</span>
|
||||
<span class="sd"> </span>
|
||||
<span class="sd"> .. math ::</span>
|
||||
<span class="sd"> x_{i+1} = x_i + \Delta x</span>
|
||||
|
||||
<span class="sd"> .. seealso::</span>
|
||||
<span class="sd"> :meth:`dxdt_Dt` for :math:`\Delta x`</span>
|
||||
<span class="sd"> x_{i,j=0} = x_{i}</span>
|
||||
<span class="sd"> </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"> </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} < \varepsilon</span>
|
||||
<span class="sd"> </span>
|
||||
<span class="sd"> .. math ::</span>
|
||||
<span class="sd"> x_{i} = x_{i,j=\text{end}}</span>
|
||||
<span class="sd"> """</span>
|
||||
<span class="n">xi</span> <span class="o">=</span> <span class="n">x0</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="n">Dx</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">xi</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">xi1</span> <span class="o">=</span> <span class="n">x0</span> <span class="o">+</span> <span class="n">Dx</span> <span class="c1"># Approximate solution at next value of x</span>
|
||||
<span class="n">residuum</span> <span class="o">=</span> <span class="n">norm</span><span class="p">(</span><span class="n">xi1</span><span class="o">-</span><span class="n">xi</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">xi1</span><span class="p">)</span>
|
||||
<span class="n">xi</span> <span class="o">=</span> <span class="n">xi1</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"><</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">xi</span><span class="p">,</span> <span class="n">iterations</span></div>
|
||||
<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="i1n"><a class="viewcode-back" href="../../numerical.html#numerical.ode.i1n">[docs]</a><span class="k">def</span> <span class="nf">i1n</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="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="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">Dx</span> <span class="o">=</span> <span class="n">dxdt_Dt</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="o">+</span><span class="mi">1</span><span class="p">],</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">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="n">iterations</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">fixed_point_iteration</span><span class="p">(</span><span class="n">Dx</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">max_iterations</span><span class="p">,</span> <span class="n">tol</span><span class="p">,</span> <span class="n">verbose</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">'Numerical integration of ODE using implicite first-order method (Euler) was successful.'</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="i1"><a class="viewcode-back" href="../../numerical.html#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>
|
||||
<div class="viewcode-block" id="i1"><a class="viewcode-back" href="../../numerical.html#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">"""Implicite first-order method / backward Euler method.</span>
|
||||
|
||||
<span class="sd"> :param f: the function to solve</span>
|
||||
@@ -286,7 +293,8 @@
|
||||
<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"> :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>
|
||||
@@ -298,26 +306,24 @@
|
||||
<span class="sd"> """</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="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="c1"># x(i+1) = x(i) + f(x(i+1), t(i+1)), exact value of f(x(i+1), t(i+1)) is not</span>
|
||||
<span class="c1"># available therefor using Newton-Raphson method</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="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">xi</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">xi1</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="c1"># Approximate solution at next value of x</span>
|
||||
<span class="n">residuum</span> <span class="o">=</span> <span class="n">norm</span><span class="p">(</span><span class="n">xi1</span><span class="o">-</span><span class="n">xi</span><span class="p">)</span><span class="o">/</span><span class="n">norm</span><span class="p">(</span><span class="n">xi1</span><span class="p">)</span>
|
||||
<span class="n">xi</span> <span class="o">=</span> <span class="n">xi1</span>
|
||||
<span class="k">if</span> <span class="n">residuum</span> <span class="o"><</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">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="p">,</span>
|
||||
<span class="n">tol</span><span class="p">,</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">'Numerical integration of ODE using implicite first-order method (Euler) was successful.'</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Numerical integration of ODE using implicite '</span> <span class="o">+</span>
|
||||
<span class="s1">'first-order method (Euler) was successful.'</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="../../numerical.html#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>
|
||||
<div class="viewcode-block" id="newmark_newtonraphson"><a class="viewcode-back" href="../../numerical.html#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">"""Newmark method.</span>
|
||||
|
||||
<span class="sd"> :param f: the function to solve</span>
|
||||
@@ -330,7 +336,8 @@
|
||||
<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"> :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>
|
||||
@@ -346,9 +353,9 @@
|
||||
<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="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>
|
||||
|
||||
@@ -361,7 +368,8 @@
|
||||
<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">#x11 = x[i,:] + dxdt*Dt # Approximate solution at next value of x</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>
|
||||
@@ -370,7 +378,8 @@
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s1">'divergiert'</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">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>
|
||||
|
||||
@@ -384,11 +393,13 @@
|
||||
<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">'Numerical integration of ODE using explicite newmark method was successful.'</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Numerical integration of ODE using explicite '</span> <span class="o">+</span>
|
||||
<span class="s1">'newmark method was successful.'</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_rdk"><a class="viewcode-back" href="../../numerical.html#numerical.ode.newmark_newtonraphson_rdk">[docs]</a><span class="k">def</span> <span class="nf">newmark_newtonraphson_rdk</span><span class="p">(</span><span class="n">fnm</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">maxIterations</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>
|
||||
<div class="viewcode-block" id="newmark_newtonraphson_rdk"><a class="viewcode-back" href="../../numerical.html#numerical.ode.newmark_newtonraphson_rdk">[docs]</a><span class="k">def</span> <span class="nf">newmark_newtonraphson_rdk</span><span class="p">(</span><span class="n">fnm</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">"""Newmark method.</span>
|
||||
|
||||
<span class="sd"> :param f: the function to solve</span>
|
||||
@@ -401,7 +412,8 @@
|
||||
<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"> :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>
|
||||
@@ -417,13 +429,14 @@
|
||||
<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="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">fnm</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">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">fnm</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>
|
||||
@@ -432,13 +445,16 @@
|
||||
<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">maxIterations</span><span class="p">):</span> <span class="c1"># Fixed-point iteration</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">#x11 = x[i,:] + dxdt*Dt # Approximate solution at next value of x</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">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>
|
||||
@@ -454,7 +470,8 @@
|
||||
<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">'Numerical integration of ODE using explicite newmark method was successful.'</span><span class="p">)</span>
|
||||
<span class="nb">print</span><span class="p">(</span><span class="s1">'Numerical integration of ODE using explicite '</span> <span class="o">+</span>
|
||||
<span class="s1">'newmark method was successful.'</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>
|
||||
|
||||
Reference in New Issue
Block a user