Ho trovato finalmente il libro che spiega il metodo di newton per approssimare le radici di una funzione.
codice:
double Newton(Function &f, const double epsilon)
{
double x0, y;
double fx0, D1fx0;
x0 = f.APointInDomain ();
do{
fx0 = f.apply(x0);
D1fx0 = (f.apply(x0+epsilon) - fx0) / epsilon;
//la precisione del double non è infinita...
if(D1fx0 == 0.0)D1fx0=epsilon;
y = x0 - (fx0 / D1fx0);
x0 = y;
}while(fabs(f.apply(y))>epsilon);
return y;
}
ora io mi sono scritto una classe Function per rendere generica la cosa, (come se lavorassi su un linguaggio di programmazione funzionale), se però ti interessa solo la radice cubica, sostituisci f.apply(x0) con x*x*x-x0 , ovvero: x al cubo meno x con zero, che è il polinomio di cui devi trovare la radice (essendo dispari ha solo 1 radice)
L'algoritmo è il seguente:
data una funzione f ed un valore e (la tolleranza di approssimazione) restituisce una radice della funzione.
siano f' la derivata prima di f, x0 un punto, r la radice della funzione
l'algoritmo inizia scegliendo un x0 nel dominio della funzione, dove la derivata sia diversa da zero, nel caso della radice cubica di x puoi scegliere x0=x
calcolo f(x0) e f'(x0)
il primo risultato sarà y = x0 - f(x0) / f'(x0)
se f(y) = 0 allora y = r
se |f(y)| < e allora y differisce da r per meno di e
altrimenti x0 = y e riinizio
piu epsilon è piccolo maggiore sarà la precisione del risultato, ma sarà maggiore anche il tempo di elaborazione, ed inoltre se epsilon è troppo piccolo potrebbe addirittura capitare che la derivata calcolata risulti zero a causa della precisione finita del tipo double nei computer.
in pratica epsilon=1E-10 dovrebbe andare bene nella maggior parte dei casi (hai un risultato preciso entro la 10a cifra decimale)
codice:
#include <iostream>
#define abs(x) ((x)<0.0?(-(x)):(x))
using namespace std;
double sqr3(const double x)
{
double x0, y;
double fx0, D1fx0;
if(x==0.0)return 0.0;
x0 = x;
do{
fx0 = x0 * x0 * x0 - x;
D1fx0 = 3 * x0 * x0;
y = x0 - (fx0 / D1fx0);
x0 = y;
}while(abs(y*y*y-x)>1e-10);
return y;
}
int main()
{
cout << sqr3(-8.0);
return 0;
}