
import numpy as np
import scipy


from scipy import integrate
help(integrate.quad)


from math import cos, exp, pi
from scipy.integrate import quad

# funzione da integrare
def f(x):
    return exp(cos(-2 * x * pi))

# chiamata a quad
res, err = quad(f, -2, 2)

print(f"The numerical result is {res:f} (+-{err:g})")


import numpy as np
def f_N(x):
    return np.exp(np.cos(-2 * x * np.pi))

res, err = quad(f_N, -2, 2)
print(f"The numerical result is {res:f} (+-{err:g})")


def my_f(x):
    return 1/x/x

res, err = quad(my_f, 1, np.inf)
print(f"The numerical result is {res:f} (+-{err:g})")


# https://www.math.ubc.ca/~pwalls/math-python/integration/trapezoid-rule/

import numpy as np

def my_trapz(f,a,b,N=50):
    '''Approximate the integral of f(x) from a to b by the trapezoid rule.

    The trapezoid rule approximates the integral \int_a^b f(x) dx by the sum:
    (dx/2) \sum_{k=1}^N (f(x_k) + f(x_{k-1}))
    where x_k = a + k*dx and dx = (b - a)/N.

    Parameters
    ----------
    f : function
        Vectorized function of a single variable
    a , b : numbers
        Interval of integration [a,b]
    N : integer
        Number of subintervals of [a,b]

    Returns
    -------
    float
        Approximation of the integral of f(x) from a to b using the
        trapezoid rule with N subintervals of equal length.

    Examples
    --------
    >>> trapz(np.sin,0,np.pi/2,1000)
    0.9999997943832332
    '''
    x = np.linspace(a,b,N+1) # N+1 points make N subintervals
    y = f(x)
    y_right = y[1:] # right endpoints
    y_left = y[:-1] # left endpoints
    dx = (b - a)/N
    T = (dx/2) * np.sum(y_right + y_left)
    return T


def vf(x):
    return np.exp(np.cos(-2 * x * np.pi))

res1 = my_trapz(vf,-2,2)
print(f"The numerical result is {res1:f}")


from scipy.integrate import trapz
#help(trapz)


x = np.linspace(-2,2,101)
y = vf(x)
res2 = trapz(y,x)
res2


from scipy.integrate import simps
#help(trapz)


x = np.linspace(-2,2,101)
y = vf(x)
res2 = simps(y,x)
res2


from scipy.integrate import odeint
import numpy as np

def f(y, t):
    """this is the rhs of the ODE to integrate, i.e. dy/dt=f(y,t)"""
    return -2 * y * t

y0 = 1             # initial value
a = 0              # integration limits for t
b = 2

t = np.arange(a, b, 0.01)  # values of t for
                          # which we require
                          # the solution y(t)
y = odeint(f, y0, t)  # actual computation of y(t)

import matplotlib.pyplot as plt         # plotting of results
fig, ax = plt.subplots()
ax.plot(t,y)
ax.set_xlabel('t')
ax.set_ylabel('y(t)')


# https://www.math.ubc.ca/~pwalls/math-python/differential-equations/first-order/

def odeEuler(f,y0,t):
    '''Approximate the solution of y'=f(y,t) by Euler's method.

    Parameters
    ----------
    f : function
        Right-hand side of the differential equation y'=f(t,y), y(t_0)=y_0
    y0 : number
        Initial value y(t0)=y0 wher t0 is the entry at index 0 in the array t
    t : array
        1D NumPy array of t values where we approximate y values. Time step
        at each iteration is given by t[n+1] - t[n].

    Returns
    -------
    y : 1D NumPy array
        Approximation y[n] of the solution y(t_n) computed by Euler's method.
    '''
    y = np.zeros(len(t))
    y[0] = y0
    for n in range(0,len(t)-1):
        y[n+1] = y[n] + f(y[n],t[n])*(t[n+1] - t[n])
    return y


y_E = odeEuler(f,y0,t)


import matplotlib.pyplot as plt         # plotting of results
fig, ax = plt.subplots()
ax.plot(t,y,label='sol by odeint')
ax.plot(t,y_E,c='r',label='sol by odeEuler')
ax.set_xlabel('t')
ax.set_ylabel('y(t)')
ax.legend()


from scipy.optimize import bisect

def f(x):
    """returns f(x)=x^3-2x^2. Has roots at
    x=0 (double root) and x=2"""
    return x ** 3 - 2 * x ** 2

# main program starts here
x = bisect(f, 1.5, 3, xtol=1e-6)

print(f"The root x is approximately x={x:14.12g},\n"
      f"the error is less than 1e-6.")
print(f"The exact error is {2 - x:g}.")


from scipy.optimize import fsolve

def f(x):
    return x ** 3 - 2 * x ** 2

xstart = 3
x = fsolve(f, xstart)      # one root is at x=2.0
                           # search starts at x=3
                           # fsolve returns a numpy array

print(f"The root x is approximately x={x[0]:21.19g}")
print(f"The exact error is {2 - x[0]:g}.")


import numpy as np
from scipy.optimize import curve_fit


def f(x, a, b, c):
    """Fit function y=f(x,p) with parameters p=(a,b,c). """
    return a * np.exp(- b * x) + c

#create fake data
x = np.linspace(0, 4, 50)
y = f(x, a=2.5, b=1.3, c=0.5)
#add noise
yi = y + 0.2 * np.random.normal(size=len(x))

#call curve fit function
popt, pcov = curve_fit(f, x, yi)
# extract fit parameters
af, bf, cf = popt
print(f"Optimal parameters are af={af:g}, bf={bf:g}, and cf={cf:g}.")

# best fit curve
yfitted = f(x,af,bf,cf)

#plot
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6,4))
ax.scatter(x, yi, marker='o', label='data $y_i$')
ax.plot(x, yfitted, c='r', label='fit $f(x_i)$')
ax.set_xlabel('x')
ax.legend()


import numpy as np
from scipy.optimize import fmin
import matplotlib.pyplot as plt

def f(x):
    return np.cos(x) - 3 * np.exp( -(x - 0.2) ** 2)

# find minima of f(x),
# starting from 1.0 and 2.0 respectively
minimum1 = fmin(f, 1.0)
print("Start search at x=1., minimum is", minimum1)
minimum2 = fmin(f, 2.0)
print("Start search at x=2., minimum is", minimum2)

# plot function
x = np.arange(-10, 10, 0.1)
y = f(x)

fig, ax = plt.subplots(figsize=(6,4))
ax.set_xlabel('x')
ax.grid(b=True,which='both')
ax.set_xlim(-5.,5.)
ax.set_ylim(-2.2,0.5)

ax.plot(x, y, label='$\cos(x)-3e^{-(x-0.2)^2}$')
# add minimum1 to plot
ax.scatter(minimum1, f(minimum1), marker='v', c='r', label='minimum 1')
# add start1 to plot
ax.scatter(1.0, f(1.0), marker='o', c='r', label='start 1')

# add minimum2 to plot
ax.scatter(minimum2,f(minimum2), marker='v', c='g', label='minimum 2')
# add start2 to plot
ax.scatter(2.0,f(2.0), marker='o', c='g',label='start 2')

ax.legend(loc='lower left')

Metodi Numerici in Python (scipy)¶

SciPy¶

Integrazione numerica¶

Integrazione su un intervallo infinito¶

Integrazione con metodo dei trapezi¶

Integrazione con metodo dei trapezi usando scipy.integrate.trapz¶

Integrazione con metodo di Simpson usando scipy.integrate.simps¶

Esercizio: integrare una funzione¶

Esercizio: fare il grafico prima di integrare¶

Risolvere equazioni differenziali ordinarie¶

Soluzione con il metodo di Eulero¶

Esercizio: usare `odeint`¶

Ricerca delle radici¶

Ricerca delle radici usando il metodo della bisezione¶

Esercizio: trovare le radici usando il metodo di bisezione¶

Ricerca delle radici usando la funzione `fsolve`¶

Fit di un set di dati¶

Ottimizzazione¶

Valori ritornati da `fmin`¶

Altri metodi numerici¶

Metodi Numerici in Python (scipy)¶

SciPy¶

Integrazione numerica¶

Integrazione su un intervallo infinito¶

Integrazione con metodo dei trapezi¶

Integrazione con metodo dei trapezi usando scipy.integrate.trapz¶

Integrazione con metodo di Simpson usando scipy.integrate.simps¶

Esercizio: integrare una funzione¶

Esercizio: fare il grafico prima di integrare¶

Risolvere equazioni differenziali ordinarie¶

Soluzione con il metodo di Eulero¶

Esercizio: usare odeint¶

Ricerca delle radici¶

Ricerca delle radici usando il metodo della bisezione¶

Esercizio: trovare le radici usando il metodo di bisezione¶

Ricerca delle radici usando la funzione fsolve¶

Fit di un set di dati¶

Ottimizzazione¶

Valori ritornati da fmin¶

Altri metodi numerici¶

Esercizio: usare `odeint`¶

Ricerca delle radici usando la funzione `fsolve`¶

Valori ritornati da `fmin`¶