{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

10: Metodi Numerici in Python (SciPy)

\n", "\n", "\n", "
\n", "\n", "
\n", "\n", "
\n", "\n", " \n", "## 10.1 Nozioni fondamentali\n", "\n", "La libreria scipy (SCIentific PYthon) fornisce un gran numero di algoritmi numerici. In questo notebook vengono presentati quelli di uso più comune che servono a calcolare integrali, trovare zeri, massimi e minimi di funzioni, risolvere equazioni differenziali e interpolare dati.
\n", "`Scipy` è basata su `numpy`. Nella [documentazione di scipy](https://docs.scipy.org/doc/scipy/reference/tutorial/general.html) viene raccomandato di importare le due librerie separatamente." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La libreria `scipy` fornisce informazioni sulla propria struttura atrraverso il comando `help`:\n", "\n", "```python\n", "help(scipy)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'output è molto lungo. Ne mostriamo solo una parte:\n", "\n", " stats --- Statistical Functions [*]\n", " sparse --- Sparse matrix [*]\n", " lib --- Python wrappers to external libraries [*]\n", " linalg --- Linear algebra routines [*]\n", " signal --- Signal Processing Tools [*]\n", " misc --- Various utilities that don't have another home.\n", " interpolate --- Interpolation Tools [*]\n", " optimize --- Optimization Tools [*]\n", " cluster --- Vector Quantization / Kmeans [*]\n", " fftpack --- Discrete Fourier Transform algorithms [*]\n", " io --- Data input and output [*]\n", " integrate --- Integration routines [*]\n", " lib.lapack --- Wrappers to LAPACK library [*]\n", " special --- Special Functions [*]\n", " lib.blas --- Wrappers to BLAS library [*]\n", " [*] - using a package requires explicit import (see pkgload)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Per trovare un algoritmo per integrare una funzione, si può esplorare la libreria `integrate`:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```python\n", "import scipy.integrate\n", "\n", "scipy.integrate?\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "produces:\n", "\n", "```text\n", "=============================================\n", "Integration and ODEs (:mod:`scipy.integrate`)\n", "=============================================\n", "\n", ".. currentmodule:: scipy.integrate\n", "\n", "Integrating functions, given function object\n", "============================================\n", "\n", ".. autosummary::\n", " :toctree: generated/\n", "\n", " quad -- General purpose integration\n", " dblquad -- General purpose double integration\n", " tplquad -- General purpose triple integration\n", " nquad -- General purpose n-dimensional integration\n", " fixed_quad -- Integrate func(x) using Gaussian quadrature of order n\n", " quadrature -- Integrate with given tolerance using Gaussian quadrature\n", " romberg -- Integrate func using Romberg integration\n", " quad_explain -- Print information for use of quad\n", " newton_cotes -- Weights and error coefficient for Newton-Cotes integration\n", " IntegrationWarning -- Warning on issues during integration\n", "\n", "Integrating functions, given fixed samples\n", "==========================================\n", "\n", ".. autosummary::\n", " :toctree: generated/\n", "\n", " trapz -- Use trapezoidal rule to compute integral.\n", " cumtrapz -- Use trapezoidal rule to cumulatively compute integral.\n", " simps -- Use Simpson's rule to compute integral from samples.\n", " romb -- Use Romberg Integration to compute integral from\n", " -- (2**k + 1) evenly-spaced samples.\n", "\n", ".. seealso::\n", "\n", " :mod:`scipy.special` for orthogonal polynomials (special) for Gaussian\n", " quadrature roots and weights for other weighting factors and regions.\n", "\n", "Integrators of ODE systems\n", "==========================\n", "\n", ".. autosummary::\n", " :toctree: generated/\n", "\n", " odeint -- General integration of ordinary differential equations.\n", " ode -- Integrate ODE using VODE and ZVODE routines.\n", " complex_ode -- Convert a complex-valued ODE to real-valued and integrate.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Attenzione!

\n", " Qualunque manipolazione di una funzione nota deve partire dal grafico che guiderà i passi successivi. Nessun metodo per quanto potente è in grado di prendere sempre le decisioni corrette per una qualsiasi funzione il cui campo di esistenza e range di valori sono in linea di principio l'intero asse reale.
\n", "
    \n", "
  1. Quando si deve fare un integrale è opportuno fare il grafico della funzione integranda per controllarne il comportamento. Singolarità (per esempio valori di $x$ in cui $f(x)$ tende a più o meno infinito) oppure altri comportamenti irregolari, come quello di $f(x)=\\sin(\\frac{1}{x}$) nell'intorno di $x = 0$, sono estremamente difficili e talvolta semplicemente impossibili da trattare numericamente.\n", "
  2. La ricerca di zeri oppure di massimi e minimi richiedono un punto oppure un intervallo iniziale. Il modo più semplice per determinare questi valori è esaminare graficamente la funzione.\n", "
\n", "Confrontare la risposta di un algoritmo con un grafico è spesso, quando possibile, il modo più semplice per verificarne la correttezza o perlomeno la non manifesta scorrettezza.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 10.2 Integrazione numerica\n", "\n", "
\n", " \n", "\n", "  \n", " (Svein Linge, Hans Petter Langtangen - Programming for Computations)\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Scientific Python fornisce diverse routines di integrazione. Uno strumento di uso generale per calcolare integrali *I* del tipo\n", "\n", "$$I=\\int_a^b f(x) \\mathrm{d} x$$\n", "\n", "è la funzione `quad()` del modulo `scipy.integrate`.\n", "\n", "Prende in input la funzione *f(x)* da integrare (l'“integrando”) e gli estremi inferiore e superiore *a* and *b*. \n", "Restituisce due valori, (in una tuple): il primo è il risultato dell'integrale mentre il secondo è una stima dell'errore numerico del risultato." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import integrate\n", "help(integrate.quad)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ecco un esempio:\n", "$$I=\\int_{-2}^{2} e^{\\cos(-2 \\pi x)} \\,\\mathrm{d} x$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from math import cos, exp, pi\n", "from scipy.integrate import quad\n", "\n", "# funzione da integrare\n", "def f(x):\n", " return exp(cos(-2 * x * pi))\n", "\n", "# chiamata a quad\n", "res, err = quad(f, -2, 2)\n", "\n", "print(f\"The numerical result is {res:f} (+-{err:g})\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Usando Numpy:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "def f_N(x):\n", " return np.exp(np.cos(-2 * x * np.pi))\n", "\n", "res, err = quad(f_N, -2, 2)\n", "print(f\"The numerical result is {res:f} (+-{err:g})\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si noti che `quad()` può prendere come parametri opzionali `epsabs` e `epsrel` per aumentare o diminuire l'accuratezza del calcolo (Usate `help(quad)` per maggiori informazioni). I valori di default sono `epsabs=1.5e-8` and `epsrel=1.5e-8`. Per il prossimo esercizio, i valori di default sono sufficienti." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integrazione su un intervallo infinito" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def my_f(x):\n", " return 1/x/x\n", "\n", "res, err = quad(my_f, 1, np.inf)\n", "print(f\"The numerical result is {res:f} (+-{err:g})\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integrazione di una funzione che dipende da parametri\n", "\n", "Talvolta la funzione da integrare dipende da dei parametri. Per esempio la velocità di un corpo che si muove in un campo gravitazionale costante è dato da $v(t) = g t + v_0 $\n", "dove $v_0$ è la velocità a $t = 0$. Lo spostamento $y(t)$ può essere calcolato integrando questa relazione fra 0 e $t$
\n", "In questo caso la variabile su cui integrare deve essere al primo posto fra le variabili della funzione e i parametri addizionali vengono passati in una ntupla `args`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import quad\n", "\n", "def vel(t,v0):\n", " g = -9.81\n", " return g*t + v0\n", "\n", "v0 = 10 # m/s\n", "tf = 3 # s\n", "\n", "y,err = quad(vel,0,tf,args=(v0,))\n", "print(f\"La posizione di un corpo che parte dall'origine con v0={v0} dopo {tf} secondi è {y} m\") " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integrazione di una singolarità \"integrabile\":" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def inv_square_root(x):\n", " return 1./np.sqrt(x)\n", "\n", "res,err = quad(inv_square_root,0,1)\n", "res,err " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integrazione di una singolarità \"non integrabile\":" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def inv(x):\n", " return 1./x\n", "\n", "res,err = quad(inv,0,1)\n", "res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integrazione con metodo dei trapezi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# https://www.math.ubc.ca/~pwalls/math-python/integration/trapezoid-rule/\n", "\n", "import numpy as np\n", "\n", "def my_trapz(f,a,b,N=50):\n", " '''Approximate the integral of f(x) from a to b by the trapezoid rule.\n", "\n", " The trapezoid rule approximates the integral \\int_a^b f(x) dx by the sum:\n", " (dx/2) \\sum_{k=1}^N (f(x_k) + f(x_{k-1}))\n", " where x_k = a + k*dx and dx = (b - a)/N.\n", "\n", " Parameters\n", " ----------\n", " f : function\n", " Vectorized function of a single variable\n", " a , b : numbers\n", " Interval of integration [a,b]\n", " N : integer\n", " Number of subintervals of [a,b]\n", "\n", " Returns\n", " -------\n", " float\n", " Approximation of the integral of f(x) from a to b using the\n", " trapezoid rule with N subintervals of equal length.\n", "\n", " Examples\n", " --------\n", " >>> trapz(np.sin,0,np.pi/2,1000)\n", " 0.9999997943832332\n", " '''\n", " x = np.linspace(a,b,N+1) # N+1 points make N subintervals\n", " y = f(x)\n", " y_right = y[1:] # right endpoints\n", " y_left = y[:-1] # left endpoints\n", " dx = (b - a)/N\n", " T = (dx/2) * np.sum(y_right + y_left)\n", " return T" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def vf(x):\n", " return np.exp(np.cos(-2 * x * np.pi))\n", "\n", "res1 = my_trapz(vf,-2,2)\n", "print(f\"The numerical result is {res1:f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integrazione con metodo dei trapezi usando scipy.integrate.trapz" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import trapz\n", "#help(trapz)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'input di `trapz` sono l'array delle coordinate `y` e quello delle coordinate `x` **in questo ordine**." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(-2,2,101)\n", "y = vf(x)\n", "res2 = trapz(y,x)\n", "res2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integrazione con metodo di Simpson usando scipy.integrate.simps\n", "\n", "Anche per `simps` l'input è costituito dall'array delle coordinate `y` e da quello delle coordinate `x` nell'ordine." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import simps\n", "#help(trapz)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(-2,2,101)\n", "y = vf(x)\n", "res2 = simps(y,x)\n", "res2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Imparare Facendo

\n", " \n", "
    \n", "
  1. Calcolate $I = \\int_{-1}^1 (x^2 -1) \\,dx$ usando quad, trapz, simps. Fate attenzione a quale tipo di informazione è necessario fornire nei vari casi. Variate il numero di punti utilizzati in trapz e simps.\n", "
  2. Usando la funzione di scipy quad, scrivete un programmma che calcola numericamente l'integrale: $I = \\int_0^1\\cos(2\\pi x) dx$.\n", "
  3. Integrare numericamente la funzione $ e^{-x/\\mu} $ in $x$ fra zero e infinito per $\\mu = 1,\\, 2,\\, 3$.\n", "
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 10.3 Risolvere equazioni differenziali ordinarie\n", "\n", "
\n", " \n", "\n", "  \n", " (Svein Linge, Hans Petter Langtangen - Programming for Computations)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si può facilmente dimostrare che una equazione differenziale ordinaria di ordine `n`, cioè che coinvolge solo derivate fino all'ordine `n` della funzione incognita, può essere riscritta come un sistema di `n` equazioni differenziali di primo grado.\n", "
 
\n", "Per risolvere equazioni differenziali ordinarie del tipo\n", "$$\\frac{\\mathrm{d}y}{\\mathrm{d}t}(t) = f(y,t)$$\n", "\n", "con condizione iniziale $y(t_0)=y_0$, oppure un sistema di equazioni tali equazioni, si può usare la funzione `odeint` di `scipy`. Ecco un esempio (`useodeint.py`) per determinare\n", "\n", "$$y(t) \\quad \\mathrm{for}\\quad t\\in[0,2]$$\n", "data l'equazione differenziale:\n", "$$\\frac{\\mathrm{d}y}{\\mathrm{d}t}(t) = -2yt \\quad \\mathrm{with} \\quad y(0)=1.$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.integrate import odeint\n", "import numpy as np\n", "\n", "def f(y, t):\n", " \"\"\"this is the rhs of the ODE to integrate, i.e. dy/dt=f(y,t)\"\"\"\n", " return -2 * y * t\n", "\n", "y0 = 1 # initial value\n", "a = 0 # integration limits for t\n", "b = 2\n", "\n", "t = np.arange(a, b, 0.01) # values of t for\n", " # which we require\n", " # the solution y(t)\n", "y = odeint(f, y0, t) # actual computation of y(t)\n", "\n", "import matplotlib.pyplot as plt # plotting of results\n", "fig, ax = plt.subplots()\n", "ax.plot(t,y)\n", "ax.set_xlabel('t')\n", "ax.set_ylabel('y(t)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Soluzione con il metodo di Eulero" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# https://www.math.ubc.ca/~pwalls/math-python/differential-equations/first-order/\n", "\n", "def odeEuler(f,y0,t):\n", " '''Approximate the solution of y'=f(y,t) by Euler's method.\n", "\n", " Parameters\n", " ----------\n", " f : function\n", " Right-hand side of the differential equation y'=f(t,y), y(t_0)=y_0\n", " y0 : number\n", " Initial value y(t0)=y0 wher t0 is the entry at index 0 in the array t\n", " t : array\n", " 1D NumPy array of t values where we approximate y values. Time step\n", " at each iteration is given by t[n+1] - t[n].\n", "\n", " Returns\n", " -------\n", " y : 1D NumPy array\n", " Approximation y[n] of the solution y(t_n) computed by Euler's method.\n", " '''\n", " y = np.zeros(len(t))\n", " y[0] = y0\n", " for n in range(0,len(t)-1):\n", " y[n+1] = y[n] + f(y[n],t[n])*(t[n+1] - t[n])\n", " return y\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y_E = odeEuler(f,y0,t)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt # plotting of results\n", "fig, ax = plt.subplots()\n", "ax.plot(t,y,label='sol by odeint')\n", "ax.plot(t,y_E,c='r',label='sol by odeEuler')\n", "ax.set_xlabel('t')\n", "ax.set_ylabel('y(t)')\n", "ax.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il comando `odeint` può prendere diversi parametri opzionali per modificare l'errore di default nell'integrazione (e per produrre output addizionale che può essere utile per il debugging). Usate il comando help per farvene un'idea:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```python\n", "help(scipy.integrate.odeint)\n", "```\n", "\n", "restituisce:\n", "\n", "```\n", "Help on function odeint in module scipy.integrate.odepack:\n", "\n", "odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)\n", " Integrate a system of ordinary differential equations.\n", " \n", " Solve a system of ordinary differential equations using lsoda from the\n", " FORTRAN library odepack.\n", " \n", " Solves the initial value problem for stiff or non-stiff systems\n", " of first order ode-s::\n", " \n", " dy/dt = func(y, t0, ...)\n", " \n", " where y can be a vector.\n", " \n", " *Note*: The first two arguments of ``func(y, t0, ...)`` are in the\n", " opposite order of the arguments in the system definition function used\n", " by the `scipy.integrate.ode` class.\n", " \n", " Parameters\n", " ----------\n", " func : callable(y, t0, ...)\n", " Computes the derivative of y at t0.\n", " y0 : array\n", " Initial condition on y (can be a vector).\n", " t : array\n", " A sequence of time points for which to solve for y. The initial\n", " value point should be the first element of this sequence.\n", " args : tuple, optional\n", " Extra arguments to pass to function.\n", " Dfun : callable(y, t0, ...)\n", " Gradient (Jacobian) of `func`.\n", " col_deriv : bool, optional\n", " True if `Dfun` defines derivatives down columns (faster),\n", " otherwise `Dfun` should define derivatives across rows.\n", " full_output : bool, optional\n", " True if to return a dictionary of optional outputs as the second output\n", " printmessg : bool, optional\n", " Whether to print the convergence message\n", " \n", " Returns\n", " -------\n", " y : array, shape (len(t), len(y0))\n", " Array containing the value of y for each desired time in t,\n", " with the initial value `y0` in the first row.\n", " infodict : dict, only returned if full_output == True\n", " Dictionary containing additional output information\n", " \n", " ======= ============================================================\n", " key meaning\n", " ======= ============================================================\n", " 'hu' vector of step sizes successfully used for each time step.\n", " 'tcur' vector with the value of t reached for each time step.\n", " (will always be at least as large as the input times).\n", " 'tolsf' vector of tolerance scale factors, greater than 1.0,\n", " computed when a request for too much accuracy was detected.\n", " 'tsw' value of t at the time of the last method switch\n", " (given for each time step)\n", " 'nst' cumulative number of time steps\n", " 'nfe' cumulative number of function evaluations for each time step\n", " 'nje' cumulative number of jacobian evaluations for each time step\n", " 'nqu' a vector of method orders for each successful step.\n", " 'imxer' index of the component of largest magnitude in the\n", " weighted local error vector (e / ewt) on an error return, -1\n", " otherwise.\n", " 'lenrw' the length of the double work array required.\n", " 'leniw' the length of integer work array required.\n", " 'mused' a vector of method indicators for each successful time step:\n", " 1: adams (nonstiff), 2: bdf (stiff)\n", " ======= ============================================================\n", " \n", " Other Parameters\n", " ----------------\n", " ml, mu : int, optional\n", " If either of these are not None or non-negative, then the\n", " Jacobian is assumed to be banded. These give the number of\n", " lower and upper non-zero diagonals in this banded matrix.\n", " For the banded case, `Dfun` should return a matrix whose\n", " rows contain the non-zero bands (starting with the lowest diagonal).\n", " Thus, the return matrix `jac` from `Dfun` should have shape\n", " ``(ml + mu + 1, len(y0))`` when ``ml >=0`` or ``mu >=0``.\n", " The data in `jac` must be stored such that ``jac[i - j + mu, j]``\n", " holds the derivative of the `i`th equation with respect to the `j`th\n", " state variable. If `col_deriv` is True, the transpose of this\n", " `jac` must be returned.\n", " rtol, atol : float, optional\n", " The input parameters `rtol` and `atol` determine the error\n", " control performed by the solver. The solver will control the\n", " vector, e, of estimated local errors in y, according to an\n", " inequality of the form ``max-norm of (e / ewt) <= 1``,\n", " where ewt is a vector of positive error weights computed as\n", " ``ewt = rtol * abs(y) + atol``.\n", " rtol and atol can be either vectors the same length as y or scalars.\n", " Defaults to 1.49012e-8.\n", " tcrit : ndarray, optional\n", " Vector of critical points (e.g. singularities) where integration\n", " care should be taken.\n", " h0 : float, (0: solver-determined), optional\n", " The step size to be attempted on the first step.\n", " hmax : float, (0: solver-determined), optional\n", " The maximum absolute step size allowed.\n", " hmin : float, (0: solver-determined), optional\n", " The minimum absolute step size allowed.\n", " ixpr : bool, optional\n", " Whether to generate extra printing at method switches.\n", " mxstep : int, (0: solver-determined), optional\n", " Maximum number of (internally defined) steps allowed for each\n", " integration point in t.\n", " mxhnil : int, (0: solver-determined), optional\n", " Maximum number of messages printed.\n", " mxordn : int, (0: solver-determined), optional\n", " Maximum order to be allowed for the non-stiff (Adams) method.\n", " mxords : int, (0: solver-determined), optional\n", " Maximum order to be allowed for the stiff (BDF) method.\n", " \n", " See Also\n", " --------\n", " ode : a more object-oriented integrator based on VODE.\n", " quad : for finding the area under a curve.\n", " \n", " Examples\n", " --------\n", " The second order differential equation for the angle `theta` of a\n", " pendulum acted on by gravity with friction can be written::\n", " \n", " theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0\n", " \n", " where `b` and `c` are positive constants, and a prime (') denotes a\n", " derivative. To solve this equation with `odeint`, we must first convert\n", " it to a system of first order equations. By defining the angular\n", " velocity ``omega(t) = theta'(t)``, we obtain the system::\n", " \n", " theta'(t) = omega(t)\n", " omega'(t) = -b*omega(t) - c*sin(theta(t))\n", " \n", " Let `y` be the vector [`theta`, `omega`]. We implement this system\n", " in python as:\n", " \n", " >>> def pend(y, t, b, c):\n", " ... theta, omega = y\n", " ... dydt = [omega, -b*omega - c*np.sin(theta)]\n", " ... return dydt\n", " ...\n", " \n", " We assume the constants are `b` = 0.25 and `c` = 5.0:\n", " \n", " >>> b = 0.25\n", " >>> c = 5.0\n", " \n", " For initial conditions, we assume the pendulum is nearly vertical\n", " with `theta(0)` = `pi` - 0.1, and it initially at rest, so\n", " `omega(0)` = 0. Then the vector of initial conditions is\n", " \n", " >>> y0 = [np.pi - 0.1, 0.0]\n", " \n", " We generate a solution 101 evenly spaced samples in the interval\n", " 0 <= `t` <= 10. So our array of times is:\n", " \n", " >>> t = np.linspace(0, 10, 101)\n", " \n", " Call `odeint` to generate the solution. To pass the parameters\n", " `b` and `c` to `pend`, we give them to `odeint` using the `args`\n", " argument.\n", " \n", " >>> from scipy.integrate import odeint\n", " >>> sol = odeint(pend, y0, t, args=(b, c))\n", " \n", " The solution is an array with shape (101, 2). The first column\n", " is `theta(t)`, and the second is `omega(t)`. The following code\n", " plots both components.\n", " \n", " >>> import matplotlib.pyplot as plt\n", " >>> plt.plot(t, sol[:, 0], 'b', label='theta(t)')\n", " >>> plt.plot(t, sol[:, 1], 'g', label='omega(t)')\n", " >>> plt.legend(loc='best')\n", " >>> plt.xlabel('t')\n", " >>> plt.grid()\n", " >>> plt.show()\n", "\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Imparare Facendo

\n", " \n", "Scrivete un programma che calcola la soluzione *y(t)* della ODE seguente, usando l'algoritmo `odeint`,\n", " $$\\frac{\\mathrm{d}y}{\\mathrm{d}t} = -\\exp(-t)(10\\sin(10t)+\\cos(10t))$$\n", " da $t=0$ a $t = 10$. Il valore iniziale è $y(0)=1$.
\n", " Mostrate graficamente la soluzione, valutandola nei punti $t=0$, $t=0.01$, $t=0.02$, ..., $t=9.99$, $t=10$.
\n", " Aiutino: una parte della soluzione $y(t)$ è presentata nella figura seguente.\n", "
\n", "\n", "
\n", " \"image\"\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## 10.4 Ricerca delle radici\n", "\n", "Cercare una $x$ tale che $f(x)=0$\n", "si chiama *ricerca delle radici* di $f$. Si noti che un problema del tipo $g(x)=h(x)$ può essere riformulato come $f(x)=g(x)−h(x)=0$.\n", "\n", "Il modulo `optimize` di `scipy` fornisce diversi strumenti per la ricerca delle radici.\n", "\n", "


\n", "\n", "
\n", " \n", "\n", "  \n", " (Svein Linge, Hans Petter Langtangen - Programming for Computations)\n", "\n", "
\n", " \n", "### Ricerca delle radici usando il metodo della bisezione \n", "\n", "L'algoritmo `bisect` è (i) robusto e (ii) concettualmente molto semplice (ma lento).\n", "\n", "Supponiamo di dover calcolare le radici di $f(x)= x^3 − 2 x^2$. Questa funzione ha una radice (doppia) in $x = 0$ e un'altra fra $x = 1.5$ (dove $f(1.5) = − 1.125$) e $x = 3$ (dove $f(3) = 9$). È facile vedere che questa altra radice si trova in $x = 2$. Ecco il programma che determina la radice numericamente:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import bisect\n", "\n", "def f(x):\n", " \"\"\"returns f(x)=x^3-2x^2. Has roots at\n", " x=0 (double root) and x=2\"\"\"\n", " return x ** 3 - 2 * x ** 2\n", "\n", "# main program starts here. Typically, the range in which to perform the search is determined from a plot.\n", "x = bisect(f, 1.5, 3, xtol=1e-6)\n", "\n", "print(f\"The root x is approximately x={x:14.12g},\\n\"\n", " f\"the error is less than 1e-6.\")\n", "print(f\"The exact error is {2 - x:g}.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il metodo `bisect()` richiede obbligatoriamente tre argomenti: (i) la funzione *f(x)*, (ii) il limite inferiore *a* (che abbiamo scelto uguale a 1.5 nell'esempio) e (ii) il limite superiore *b* (scelto uguale a 3) dell'intervallo in cui eseguire la procedura. Il parametro opzionale `xtol` determina l'errore massimo del metodo.\n", "\n", "Uno dei presupposti del metodo di bisezione è che l'intervallo \\[*a*, *b*\\] sia scelto in modo tale che la funzione abbia in *a* segno opposto a quello che ha in *b* in modo che nell'intervallo, se la funzione è continua, cada almeno una radice." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Imparare Facendo

\n", " \n", "Scrivete un programma chiamato `sqrttwo.py` per cercare la radice *x* della funzione $f(x)=2 − x^2$ usando il metodo di bisezione. Scegliete come tolleranza per l'approssimazione alla radice di 10−8.\n", "
    \n", "
  1. Giustificate la scelta dell'intervallo iniziale $[a, b]$ per la ricerca: che valori evte scelto per *a* e *b*? Perchè?\n", "
  2. Esaminate i risultati:\n", "
      \n", "
    1. Che valore restituisce l'algoritmo di bisezione per la radice *x*?\n", "
    2. Calcolate il valore di $\\sqrt{2}$ usando math.sqrt(2) e confrontatelo il risultato precedente. Quanto è grande l'errore assoluto? Come si confronta con xtol?\n", "
    \n", "
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Ricerca delle radici usando la funzione `fsolve`\n", "\n", "Un algoritmo per la ricerca delle radici che è (spesso) migliore (nel senso di “più efficiente”) di quello di bisezione è codificato nella funzione `fsolve()` che funziona anche per problemi a più dimensioni. Questo algoritmo richiede solamente un punto di partenza vicino a dove ci si aspetta che ci sia una radice (Non è però detto che il metodo converga).\n", "\n", "Ecco un esempio:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import fsolve\n", "\n", "def f(x):\n", " return x ** 3 - 2 * x ** 2\n", "\n", "xstart = 3\n", "x = fsolve(f, xstart) # one root is at x=2.0\n", " # search starts at x=3\n", " # fsolve returns a numpy array\n", "\n", "print(f\"The root x is approximately x={x[0]:21.19g}\")\n", "print(f\"The exact error is {2 - x[0]:g}.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il valore restituito da `fsolve` è un array di numpy di lunghezza $n$ per un problema di ricerca di radici con $n$ variabili. Nell'esempio precedente, $n = 1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Imparare Facendo

\n", " \n", "
    \n", "
  1. Studiate il grafico della funzione y(x) = np.cos(2*x) - np.cosh(-x+1) + 2, cambiando l'intervallo di variazione dell'ascissa fino determinare il numero e la posizione approssimata degli zeri. Trovate il valore degli zeri usando fsolve \n", "
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 10.5 Fit di un set di dati\n", "\n", "Scipy fornisce una funzione piuttosto flessibile (basato sull'algoritmo di Levenburg-Marquardt), `scipy.optimize.curve_fit`, per interpolare un set di dati. L'assunzione è che vengano dati un set di punti\n", "$x_1, x_2,\\cdots,x_N$, i corrispondenti valori $y_i$ e una dipendenza funzionale $y=f(x,\\vec{p})$.\n", "Per fare un esempio, il numero di atomi non decaduti in un campione di una sostanza radioattiva segue la legge\n", "$$ N(t,N_0,\\tau) = N_0\\,\\exp\\left(-\\frac{t}{\\tau}\\right).$$\n", "Si vuole determinare i parametri $\\vec{p}=(p_1, p_2, \\ldots,p_k)$ in modo che $r$, la somma degli scarti quadratici fra la curva e i dati, sia la più piccola possibile:\n", "\n", "$$r = \\sum\\limits_{i=1}^N \\left(y_i - f(x_i, \\vec{p})\\right)^2$$\n", "\n", "Questo approccio è particolarmente utile quando i dati sono affetti da *rumore*: per ogni coppia $x_i,y_i$ è presente un termine di errore (ignoto) $\\epsilon_i$, tale che $y_i=f(x_i,\\vec{p})+\\epsilon_i$.\n", "\n", "Un esempio per chiarire: assumiamo di avere dei dati che sappiamo essere descritti dalla funzione:\n", "$$f(x,\\vec{p}) = a \\exp(-b x) + c,$$\n", "che dipende dai parametri $\\vec{p}=(a,b,c)$, che devono essere determinati usando i dati." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.optimize import curve_fit\n", "\n", "\n", "def f(x, a, b, c):\n", " \"\"\"Fit function y=f(x,p) with parameters p=(a,b,c). \"\"\"\n", " return a * np.exp(- b * x) + c\n", "\n", "#create fake data\n", "x = np.linspace(0, 4, 50)\n", "y = f(x, a=2.5, b=1.3, c=0.5)\n", "#add noise\n", "yi = y + 0.2 * np.random.normal(size=len(x))\n", "\n", "#call curve fit function\n", "popt, pcov = curve_fit(f, x, yi)\n", "# extract fit parameters\n", "af, bf, cf = popt\n", "print(f\"Optimal parameters are af={af:g}, bf={bf:g}, and cf={cf:g}.\")\n", "\n", "# best fit curve\n", "yfitted = f(x,af,bf,cf)\n", "\n", "#plot\n", "import matplotlib.pyplot as plt\n", "fig, ax = plt.subplots(figsize=(6,4))\n", "ax.scatter(x, yi, marker='o', label='data $y_i$')\n", "ax.plot(x, yfitted, c='r', label='fit $f(x_i)$')\n", "ax.set_xlabel('x')\n", "ax.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si noti che nell'esempio precedente abbiamo definito la funzione da utilizzare per il fit $y = f(x)$ usando del codice Python. Quindi possiamo utilizzare una funzione (quasi) arbitraria nel metodo `curve_fit`.\n", "\n", "La funzione `curve_fit` restituisce una tuple `popt, pcov`. Il primo elemento, `popt`, contiene l'ntupla dei parametri ottimali (OPTimal Parameters), cioè i parametri che minimizzano la somma degli scarti quadratici. Il secondo elemento contiene la matrice di covarianza di tutti i parametri (Sperimentazioni di Fisica). La diagonale fornisce la varianza della stima dei parametri. La standard deviation sulla stima dei parametri può quindi essere calcolata come\n", "```\n", "perr = np.sqrt(np.diag(pcov)).\n", "```\n", "L'algoritmo di Levenburg-Marquardt richiede, per iniziare la sua procedura, una stima iniziale dei parametri. Se questa non viene fornita, come nell'esempio precedente, il valore “1.0“ viene assunto come stima di partenza.\n", "\n", "Se l'algoritmo non riesce a fittare i dati con la funzione ipotizzata, è necessario fornire a `curve_fit` una stima migliore per i parametri iniziali. Nell'esempio precedente potremmo passare le nostre stime cambiando la linea\n", "\n", "```python\n", "popt, pcov = curve_fit(f, x, yi)\n", "```\n", "\n", "in\n", "\n", "```python\n", "popt, pcov = curve_fit(f, x, yi, p0=(2,1,0.6))\n", "```\n", "se avessimo ragione di ritenere che *a* = 2, *b* = 1 and *c* = 0.6 siano dei valori \"ragionevoli\". Una volta che la stima iniziale è \"più o meno corretta\" il fit funziona bene. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Imparare Facendo

\n", " \n", "Un set di misure ha dato il risultato\n", "\n", "array([ 1.98781717e-01, 3.24825975e-01, 2.98912925e-01, 4.93741345e-02,\n", " 3.80536122e-01, -2.05719267e-01, -7.71630572e-03, -1.51733833e-01\n", " 1.73748896e-02, 9.02509179e-03, -2.20920245e-02, 1.34501550e-02\n", " 3.06908554e-01, 2.34147443e-01, 7.75703310e-01, 4.99104499e-01\n", " 5.66648179e-01, 1.04641129e+00, 1.67507051e+00 1.77283086e+00\n", " 1.96239227e+00, 2.27163706e+00, 2.54161451e+00 3.09807009e+00\n", " 3.45291476e+00, 4.23653254e+00, 4.58924190e+00 5.08824411e+00\n", " 5.76377627e+00, 6.04351916e+00, 6.59161494e+00 7.25319841e+00\n", " 7.96917169e+00, 8.06559558e+00, 8.91595385e+00 9.78458262e+00\n", " 1.04143742e+01, 1.14280612e+01, 1.18132229e+01 1.25599940e+01\n", " 1.35895984e+01, 1.41159890e+01, 1.53848688e+01 1.63810995e+01\n", " 1.69877641e+01, 1.77683845e+01, 1.88923474e+01 1.95531287e+01\n", " 2.12809123e+01, 2.15518168e+01])\n", "\n", "ai tempi\n", "\n", "array([0. , 0.10204082, 0.20408163, 0.30612245, 0.40816327,\n", " 0.51020408, 0.6122449 , 0.71428571, 0.81632653, 0.91836735,\n", " 1.02040816, 1.12244898, 1.2244898 , 1.32653061, 1.42857143,\n", " 1.53061224, 1.63265306, 1.73469388, 1.83673469, 1.93877551,\n", " 2.04081633, 2.14285714, 2.24489796, 2.34693878, 2.44897959,\n", " 2.55102041, 2.65306122, 2.75510204, 2.85714286, 2.95918367,\n", " 3.06122449, 3.16326531, 3.26530612, 3.36734694, 3.46938776,\n", " 3.57142857, 3.67346939, 3.7755102 , 3.87755102, 3.97959184,\n", " 4.08163265, 4.18367347, 4.28571429, 4.3877551 , 4.48979592,\n", " 4.59183673, 4.69387755, 4.79591837, 4.89795918, 5. ])\n", "\n", "Trovate la parabola della forma \n", "y(t) = a*t**2 + b*t + c\n", "\n", "che meglio approssima i dati.
\n", "Fate un plot dei dati e della curva ottenuta.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 10.6 Ottimizzazione ovvero trovare massimi (e minimi)\n", "\n", "Spesso è necessario trovare il massimo o il minimo di una particolare funzione *f(x)* dove *f* è una funzione scalare mentre *x* può essere un vettore. Applicazioni tipiche sono la minimizzazione di quantità come il costo, il rischio o l'errore, oppure la massimizzazione della produttività, efficienza o profitto. Le routines di ottimizzazione tipicamente forniscono un metodo per minimizzare una funzione data: per massimizzare *f(x)* è sufficiente minimizzare *g(x)= − f(x)*.\n", "\n", "Di seguito, un esempio che mostra (i) la definizione della funzione *f(x)* da massimizzare e (ii) la chiamata a `scipy.optimize.fmin` a cui vengono passati la funzione *f* da minimizzare e un valore iniziale *x*0 da cui partire per la ricerca del minimo, e che restituisce il valore *x* per cui *f(x)* ha un minimo locale. Tipicamente, la ricerca del minimo è locale, nel senso che l'algoritmo segue il gradiente (derivata multidimensionale) nel punto in cui si trova. \n", "Facciamo come prima cosa il grafico della funzione. Il plot suggerisce che i due minimi più a destra possono essere cercati partendo dai due punti *x*0 = 1.0 e *x*0 = 2.0. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.optimize import fmin\n", "import matplotlib.pyplot as plt\n", "\n", "def f(x):\n", " return np.cos(x) - 3 * np.exp( -(x - 0.2) ** 2)\n", "\n", "\n", "# plot function \n", "x = np.arange(-10, 10, 0.1)\n", "y = f(x)\n", "\n", "fig, ax = plt.subplots(figsize=(6,4))\n", "ax.set_xlabel('x')\n", "ax.grid(visible=True,which='both')\n", "ax.set_xlim(-5.,5.)\n", "ax.set_ylim(-2.2,0.5)\n", "\n", "ax.plot(x, y, label='$\\cos(x)-3e^{-(x-0.2)^2}$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Troviamo i due minimi con due chiamate a `fmin`. Il risultato mostra come, in generale, partendo da punti diversi troviamo minimi diversi. È anche possibile esistano punti iniziali che non portano a nessun minimo. Rigeneriamo il plot della funzione, mostrando i punti iniziali delle due ricerche e i minimi ottenuti:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# find minima of f(x),\n", "# starting from 1.0 and 2.0 respectively\n", "minimum1 = fmin(f, 1.0)\n", "print(\"Start search at x=1., minimum is\", minimum1)\n", "minimum2 = fmin(f, 2.0)\n", "print(\"Start search at x=2., minimum is\", minimum2)\n", "\n", "# plot function again showing the minima that have been found\n", "x = np.arange(-10, 10, 0.1)\n", "y = f(x)\n", "\n", "fig, ax = plt.subplots(figsize=(6,4))\n", "ax.set_xlabel('x')\n", "ax.grid(visible=True,which='both')\n", "ax.set_xlim(-5.,5.)\n", "ax.set_ylim(-2.2,0.5)\n", "\n", "ax.plot(x, y, label='$\\cos(x)-3e^{-(x-0.2)^2}$')\n", "# add minimum1 to plot\n", "ax.scatter(minimum1, f(minimum1), marker='v', c='r', label='minimum 1')\n", "# add start1 to plot\n", "ax.scatter(1.0, f(1.0), marker='o', c='r', label='start 1')\n", "\n", "# add minimum2 to plot\n", "ax.scatter(minimum2,f(minimum2), marker='v', c='g', label='minimum 2')\n", "# add start2 to plot\n", "ax.scatter(2.0,f(2.0), marker='o', c='g',label='start 2')\n", "\n", "ax.legend(loc='lower left')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Come si vede, la chiamata a `fmin` produce anche dell'output addizionale che può essere utile per analizzare la procedura.\n", "\n", "##### Valori ritornati da `fmin`\n", "\n", "Si noti che la funzione `fmin` restituisce un numpy `array` che – nel caso precedente – contiene un solo numero dal momento che abbiamo una sola variabile (qui *x*) da variare. In generale, `fmin` può essre usata per trovare il minimo in uno spazio pluridimensionale. In questo caso, il numpy array contiene le coordinate del punto che minimizza la funzione obiettivo." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Imparare Facendo

\n", " \n", "
    \n", "
  1. Trovate massimi e minimi della funzione $f(x) = 2\\, x^3 + 3\\, x^2 - 12\\, x - 10$.\n", "
  2. Trovate la distanza minima fra il punto del piano (3,2) e la retta y = 3x -1.\n", "
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Altri metodi numerici\n", "-----------------------\n", "\n", "Scientific Python and Numpy forniscono molti altri algoritmi numerici: per esempio interpolazione di funzioni, trasformate di Fourier, ottimizzazione, funzioni speciali (Funzioni di Bessel etc.), generazione di numeri casuali, signal processing e filtri." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Ulteriori informazioni ed esempi:
\n", "[Scipy-lectures](http://scipy-lectures.org/index.html).
\n", "[Robert Johansson - Scientific Computing with Python](http://raw.github.com/jrjohansson/scientific-python-lectures/master/Scientific-Computing-with-Python.pdf)
\n", "[Sven Linge, Hand Petter Langtangen - Programming for Computations](https://link.springer.com/book/10.1007/978-3-319-32428-9)" ] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.12" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }