Càlcul numèric

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ó.

Deixa un comentari

L'adreça electrònica no es publicarà.

Aquest lloc utilitza Akismet per reduir els comentaris brossa. Apreneu com es processen les dades dels comentaris.