9. Equazioni algebriche non lineari#
Questa sezione si occupa della soluzione delle equazioni algebriche non lineari, distinguendo le equazioni non lineari con una sola incognita \(x \in \mathbb{R}\)
e i sistemi di equazioni non lineari con un numero di incognite pari al numero di equazione,
9.1. Equazioni non lineari#
Vengono presentati i metodi di bisezione e di Newton per la soluzione di un’equazione non lineare,
e applicati alla soluzione del problema con \(f(x) = e^x + x\), la cui derivata è nota e immediata da calcolare \(f'(x) = e^x + 1\). L’espressione della derivata verrà utilizzata nel metodo di Newton.
""" Import libaries """
import numpy as np
from time import time
"""
Example. f(x) = e^x + x
"""
# Function f and its derivative
f = lambda x: np.exp(x) + x
df = lambda x: np.exp(x) + 1
9.1.1. Metodo di bisezione#
Il metodo di bisezione per la ricerca degli zeri di una funzione continua \(F(x)\) si basa sul teorema dei valori intermedi per le funzioni continue.
Dati due numeri reali \(a\), \(b\) tali che \(f(a) \, f(b) < 0\), allora esiste un punto \(c \in (a,b)\) tale che \(f(c) = 0\).
"""
Define bisection_method_scalar() function to solve nonlinear scalar equations with bisection method
"""
def bisection_method_scalar(f, a, b, tol=1e-6, max_niter=100):
""" Function implementing the bisection method for scalar equations """
niter = 0
if ( not f(a) * f(b) < 0 ):
print("Bisection algorithm can't start, f(a)f(b)>= 0")
else:
x = .5 * (a+b)
fx = f(x)
while ( np.abs(fx) > tol and niter < max_niter ):
if ( f(x) * f(a) <= 0 ): # new range [a,c]
b = x
else: # new range [a,b]
a = x
# Update solution and residual
x = .5 * (a+b)
fx = f(x)
# Update n.iter
niter += 1
return x, np.abs(fx), niter, max_niter
""" Use bisection_method_scalar() function to solve the example """
# Find 2 values so that $f(a) f(b) < 0$
a, b = -2., 0.
t1 = time()
x, res, niter, max_niter = bisection_method_scalar(f, a, b,)
print("Bisection method summary: ")
if ( niter < max_niter ):
print(f"Convergence reached")
print(f"Sol, x = {x}")
else:
print(f"max n.iter reached without convergence")
print(f"residual : {f(x)}")
print(f"n. iterations: {niter}")
print(f"elapsed time : {time()-t1}")
Bisection method summary:
Convergence reached
Sol, x = -0.567143440246582
residual : -2.348157265297246e-07
n. iterations: 20
elapsed time : 0.0011241436004638672
9.1.2. Metodo di Newton#
Per trovare la soluzione del problema non lineare
il metodo di Newton sfrutta l’espansione in serie troncata al primo grado della funzione \(f(x)\), per scrivere
e ottenere l’incremento della soluzione \(\Delta x\) come soluzione del sistema lineare
e aggiornare la soluzione \(x^{n+1} = x^{n} + \Delta x\).
"""
Define newton_method_scalar() function to solve nonlinear scalar equations with Newton's method
"""
def newton_method_scalar(f, df, x=.0, tol=1e-6, max_niter=100):
""" Function implementing Newton's method for scalar equations """
res = f(x)
niter = 0
# Newton algorithm
while ( np.abs(res) > tol and niter < max_niter ):
# Solve linear approximation step, and update solution
dx = - res / df(x)
x += dx
#> Evaluate new residual and n. of iter
res = f(x)
niter += 1
return x, res, niter, max_niter
""" Use newton_method_scalar() function to solve the example """
# import numpy as np # already imported
# Parameters of the Newton method, for stopping criteria
# tol = 1e-6 # tolerance on the residual |f(x)| < tol
# max_niter = 10 # max n. of iterations niter > max_niter
x0 = -1.
t1 = time()
x, res, niter, max_niter = newton_method_scalar(f, df, x=x0)
print("Newton's method summary: ")
if ( niter < max_niter ):
print(f"Convergence reached")
print(f"Sol, x = {x}")
else:
print(f"max n.iter reached without convergence")
print(f"residual : {f(x)}")
print(f"n. iterations: {niter}")
print(f"elapsed time : {time()-t1}")
Newton's method summary:
Convergence reached
Sol, x = -0.567143285989123
residual : 6.927808993140161e-09
n. iterations: 3
elapsed time : 0.00038504600524902344
9.2. Sistemi di equazioni non lineari#
9.2.1. Metodo di Newton#
Il metodo di Newton sfrutta l’espansione lineare della funzione \(\mathbf{f}(\mathbf{x})\) nell’intorno di un valore \(\mathbf{x}\),
per costruire un metodo iterativo composto da due passi a ogni iterazione:
ricerca dell’incremento:
\[\mathbf{f}'(\mathbf{x}^{n}) \, \mathbf{h}^{(n)} = - \mathbf{f}(\mathbf{x}^{(n)})\]aggiornamento della soluzione
\[\mathbf{x}^{(n+1)} = \mathbf{x}^{(n)} + \mathbf{h}^{(n)} \ .\]
Esempio. Il metodo di Newton viene applicato al sistema non lineare
che può essere scritto con il formalismo matriciale come
La derivata della funzione \(\mathbf{f}(\mathbf{x})\), rispetto alla variabile indipendente \(\mathbf{x}\),
può essere rappresentata come una matrice che ha come elemento alla riga \(i\) e alla colonna \(j\) la derivata della funzione \(f_i(\mathbf{x})\) rispetto alla variabile \(x_j\), \(\left[ \mathbf{f}' \right]_{ij} = \frac{\partial f_i}{\partial x_j}\), così che l’approssimazione al primo ordine dell’incremento della funzione può essere scritto come
"""
Example. f(x) = np.array([ x[0] - x[1]
-x[0]**2 + x[1] - 1 ])
"""
# Function f and its derivative
f = lambda x: np.array([ x[0] - x[1], -x[0]**2 + x[1] + 1])
df = lambda x: np.array([[1, -1], [-2*x[0], 1]])
"""
Define newton_method_scalar() function to solve nonlinear systems of equations with Newton's method
"""
def newton_method_system(f, df, x=.0, tol=1e-6, max_niter=100):
""" Function implementing Newton's method for systems of equations """
res = f(x)
niter = 0
# Newton algorithm
while ( np.linalg.norm(res) > tol and niter < max_niter ):
# Solve linear approximation step, and update solution
dx = - np.linalg.solve(df(x), res)
x += dx
#> Evaluate new residual and n. of iter
res = f(x)
niter += 1
return x, res, niter, max_niter
""" Use newton_method_scalar() function to solve the example """
# import numpy as np # already imported
# Parameters of the Newton method, for stopping criteria
# tol = 1e-6 # tolerance on the residual |f(x)| < tol
# max_niter = 10 # max n. of iterations niter > max_niter
x0 = np.array([ -1. , 1 ])
t1 = time()
x, res, niter, max_niter = newton_method_system(f, df, x=x0)
print("Newton's method summary: ")
if ( niter < max_niter ):
print(f"Convergence reached")
print(f"Sol, x = {x}")
else:
print(f"max n.iter reached without convergence")
print(f"residual : {f(x)}")
print(f"n. iterations: {niter}")
print(f"elapsed time : {time()-t1}")
Newton's method summary:
Convergence reached
Sol, x = [-0.61803399 -0.61803399]
residual : [ 0.00000000e+00 -2.10942375e-13]
n. iterations: 4
elapsed time : 0.0007607936859130859
todo L’algoritmo di Newton trova solo una soluzione del problema. Cercare le altre soluzioni cambiando il tentativo iniziale.
todo … altro?