Aquest matí m’he entretingut a passar el càlcul (aproximat) d’un zero d’una funció del mètode de la bisecció al de Newton-Raphson.

En Python, si tenim la funció f(t) i la seva derivada f'(t) definides:

def f(t):
   ...
   return valor_funcio_a_t
def fprima(t):
   ...
   return valor_primera_derivada_a_t

l’algorisme queda tan senzill com:

t_n = 0.0
f_n = f(t_n)
while abs(f_n) > 0.0001:
    t_n = t_n - f_n / fprima(t_n)
    f_n = f(t_n)
return t_n

O bé:

t_n = 0.0
while True:
    f_n = f(t_n)
    if abs(f_n) < 0.0001:
        break
    else:
        t_n = t_n - f_n / fprima(t_n)
return t_n

O, amb versions a partir de la 3.8:

t_n = 0.0
while (f_n := f(t_n)) > 0.0001:
    t_n = t_n - f_n / fprima(t_n)
return t_n

Ho he deixat amb aquesta última versió (és la primera vegada que faig servir aquesta sintaxi de comparació i assignació combinades).

Degut a la naturalesa de la funció que avaluo no em cal prendre mesures per detectar situacions en què el mètode de càlcul no convergeix o si la derivada és zero: la meva funció és estrictament creixent. Queda breu i elegant.

Hi ha 2 paràmetres que, en general, cal adaptar a la naturalesa del problema:

• El primer valor de la variable amb què s’inicia l’algorisme. En el meu cas el resultat pot ser un número positiu o negatiu, començar amb zero és una bona tria.
• La condició de sortida (“resultat prou exacte”): la funció que avaluo pren valors de milers, un resultat en què els tres pimers decimals són zero (és a dir, la funció queda avaluada entre -0,0001 i 0,0001) és més que suficient.

També he provat el mètode de Halley (que no coneixia), però no hi guanyo pràcticament res (1 o 2 iteracions sobre 7 o 8, però a cada iteració hi ha força més feina, per calcular la segona derivada).

1/1/2022: En un cas especial he tingut derivada zero; he retocat el programa perquè tracti l’excepció.