Una versione in rete di anni fa: http://www.elegio.it/calcolatrice/trovaradici3-4.html

Questa ( 15 maggio 2013): http://www.elegio.it/calcolatrice/quarto-grado-1.html

Pagina forse utile almeno come esempio di programmazione in Javascript e in mathematichese ( http://en.wikipedia.org/wiki/Mathematica )...
Ovviamente ci sono formule migliori ma non molto usate e poco conosciute. Basta fare qualche ricerca sulla Wikipedia...

Radici di 3° e 4° grado affidabili ( abbastanza ) ...

Il calcolo numerico richiede algoritmi solidi ossia tali da dare risultati affidabili in ogni situazione, anche la meno prevedibile. L'algoritmo potrebbe venire utilizzato per milioni di volte nel corso di una qualche simulazione e... non dovrà MAI andare in crisi...

Le formule di soluzione delle equazioni di terzo e quarto grado riportate in molti testi liceali (e non) non sono abbastanza valide perché con determinati valori dei coefficienti la precisione di calcolo diventa bassa, con grande spreco di cifre significative.
Nel seguito vengono proposti algoritmi che dovrebbero favorire un uso più frequente dei polinomi di quarto grado come funzioni per l'interpolazione di funzioni complesse note solo in pochi punti. Ovviamente un polinomio di quarto grado dovrebbe consentire approssimazioni migliori di uno di secondo o addirittura di primo e il calcolo del valore della variabile indipendente in corrispondenza del quale il polinomio assume un valore prefissato non presenta, con le seguenti formule, assolutamente nessuna difficoltà o criticità numerica.
Ben altra situazione si potrebbe avere usando le formule che si trovano in certi prontuari ... probabilmente non aggiornati o usati acriticamente da SECOLI.

Noto che anche prontuari diffusi come Numerical Recipes non riportano (anno 2003, V 2.10 pag. 178 cap. 5.6) le formule per il polinomio di quarto grado disincentivando i programmatori all'uso di un' ottima approssimazione di grado elevato che pure non presenterebbe nessuna controindicazione di carattere operativo.

Noto anche che l' Handbook of Mathematical Functions di Abramowitz and Stegun a pag. 17 e pag. 18 che pure ha il merito di segnalare formule valide, non approfondisce certe questioni piuttosto cruciali per non sperperare cifre significative.

Incidentalmente vengono dati alcuni suggerimenti per l'utilizzo di Javascript come linguaggio per la didattica dell'algebra e della matematica in genere. Le formule riportate sono tutte verificabili appunto perché programmate in Javascript, popolarissimo e gratuito linguaggio interpretato. Gli algoritmi visualizzati, se i calcoli funzionano, non contengono errori di stampa poiché sono la trascrizione automatizzata delle istruzioni che vengono eseguite attivando i pulsanti.

Equazione di secondo grado con coefficienti complessi

Le funzioni usate sfruttano le capacità vettoriali del linguaggio Javascript ( alias JScript, alias EcmaScript).

Le funzioni possono fornire come risultato dei vettori e gli elementi di un vettore possono essere a loro volta vettori. Per rappresentare i numeri complessi basta dunque usare vettori di due elementi (almeno).

Nel seguito viene mostrata la funzione di calcolo di un polinomio a coefficienti complessi.

I coefficienti sono raggruppati in un vettore in cui ogni elemento è un numero complesso ossia è, a sua volta, un vettore di almeno due elementi.

.

Risolvere una equazione di secondo grado sembrerebbe questione di matematica liceale ma le insidie di carattere numerico rendono utile la disponibilità di una versione ben programmata di questo algoritmo.

Il nocciolo della questione è il fatto che la formula normalmente usata comporta l'effettuazione di due operazioni di somma: il determinante viene sommato al termine di primo grado per trovare la prima delle radici e sottratto per trovare la seconda. La somma di due quantità entrambe positive o entrambe negative non presenta difficoltà numeriche ma la differenza può provocare errori numerici anche fortissimi.

Si immagini di dover sottrarre due numeri molto vicini tra di loro (ad esempio 1.23457834 e 1.23426603). Ciascuno di essi è affetto da un errore di pochi milionesimi ( le cifre in rosso supponiamo che siano quelle errate ) ma dato che i due numeri differiscono tra loro di pochi decimillesimi l'errore della differenza sarà dell'ordine del percento !!!! ( ovvero, nel nostro caso si ha: 3.1231e-5 )

Fortunatamente la precisione non va persa se si cerca la radice di modulo maggiore e si ottiene la radice di modulo minore dividendo il termine noto (che è il prodotto delle due radici) per la radice di modulo maggiore.

Questo algoritmo è stato pensato per la ricerca delle radici di polinomi di quarto grado, che possono richiedere soluzioni di polinomi di secondo grado a coefficienti complessi; tuttavia merita di essere considerato per il suo valore didattico sia in ambito di calcolo numerico che di programmazione in Javascript (alias JScript, alias EcmaScript) per fini matematici.

Il linguaggio Javascript serve principalmente per dinamicizzare le pagine HTML ma si presta alla didattica della matematica perché è non solo gratuito ed enormemente diffuso ma anche notevolmente e sorprendentemente potente dal punto di vista del calcolo numerico.

Ovviamente non va applicato nel calcolo intensivo, dato che si tratta pur sempre di un linguaggio interpretato, ma permette una buona sinteticità nella enunciazione di algoritmi e kernel di interesse matematico.

La gestione dei numeri complessi è un esempio a sostegno di queste affermazioni.

La funzione che risolve l'equazione di secondo grado a coefficienti complessi è la seguente:

.

x2 + a x + b = 0
ar=; ai=; br=; bi=;

...

Si noti l'uso della istruzione with(...){...} che serve per scrivere formule con adeguata sinteticità.

Il linguaggio Javascript è orientato agli oggetti, come lo sono tutti i linguaggi moderni di programmazione. Nella scrittura di una formula la specifica dell'appartenenza di una proprietà o di un metodo ad un dato oggetto disturba perché la consuetudine matematica e la leggibilità vogliono che le grandezze usate nella formula siano indicate con simboli di un unico carattere, non con molti caratteri come risulta necessario quando si specifica anche l'oggetto a cui la variabile appartiene.

In Javascript le funzioni matematiche appartengono ad un oggetto predefinito di nome Math per cui volendo usare una funzione quale il valore assoluto, la radice quadrata, le funzioni trigonometriche etc.... il programmatore dovrebbe specificare, per ciascuna di esse, il fatto che sono tutte metodi dell' oggetto Math. Dunque scrivere Math.sin(0.3) oppure Math.exp(1.23) oppure Math.sqrt(0.4) e così via.
L'istruzione with(...){...} consente di sottintendere l'oggetto che... è ovviamente usato e dunque si potrà usare la normale notazione matematica avendo specificato with(Math){...}.

Si noti che la with può essere nidificata e che, in caso di conflitto, ossia se una data proprietà è definita in due diversi tipi di oggetto, prevale l'oggetto indicato con la with più interna.

Ad esempio supponiamo di definire un oggetto che abbia come sua proprietà la componente PI. Nell'oggetto Math è già definito Math.PI che è la costante π ossia vale 3.141592653589793....

Supponiamo dunque di aver scritto la seguente function che definisce un nostro oggetto di tipo ambiguo nel nome e nei fatti: 

.

Se ora eseguiamo questa funzione in cui si fa uso di un oggetto di tipo ambiguo: 

.

abbiamo la sorpresa di vedere che PI vale.... boh...

Molto rapidamente si è mostrato qui il metodo da usare per definire metodi in un oggetto non predefinito. Detto in parole povere, se si usa una funzione senza specificare la lista, anche vuota, dei suoi argomenti, in Javascript si intende che la variabile a sinistra del segno uguale non riceve il valore calcolato dalla funzione ma diventa un sinonimo della funzione ossia riceve il testo di definizione della funzione.

Per definire un oggetto di tipo personale basta scrivere una funzione che restituisca this dove questa parola chiave indica gli elementi stessi dell'oggetto che verrà indicato col nome della funzione stessa. Dunque in questo esempio si è creato un tipo di dato che ha per componenti la proprietà PI ed un metodo di nome radice che usa, ed è, un alias del metodo Math.sqrt.

Nella function provadubbio viene dichiarata la variabile insidia che è definita di tipo ambiguo (notare l'istruzione new ambiguo()). Infine con la istrizione with(insidia) abbiamo fatto in modo che l'appartenenza di radice ad insidia venisse sottintesa e pertanto, scrivendo radice(2), abbiamo attivato insidia.radice(2) ed abbiamo ottenuto 1.41421356...

Equazione di terzo grado a coefficienti reali

Anche la soluzione dell'equazione di terzo grado presenta insidie numeriche che rendono consigliabile l'uso di un algoritmo ben sperimentato dal punto di vista della minimizzazione degli errori di troncamento. Riporto qui l'algoritmo base pubblicato da Tartaglia e Cardano nel 1545 ( litigarono sulla paternità) :

http://it.wikipedia.org/wiki/Equazione_di_terzo_grado

ed una versione che risulta citata dalla wikipedia inglese:

http://en.wikipedia.org/wiki/Cubic_function 

.

x3 + p x2 + q x + r = 0
p=; q=; r=;

...

Usando invece la funzione cubicabis che usa funzioni esponenziali e trigonometriche...

...

Riporto le espressioni in stile matematico ( mi prendo la libertà di indicare il classico if(test){espressione} usando una parentesi sottolineata collegata alla condizione ossia (condizione)(espressione)  ):

a = ( p2 − 3 q ) / 9
b = ( 9 p q − p3 − 27 r ) / 54
d = b2 − a3
( d ≥ 0 ) (  
 β = ( d1/2 + | b | )1/3 b / | b |
 γ = a / β
 x1 = β + γ − p / 3
 x2 = − ( β + γ ) / 2 − p / 3 + 31/2 ( β − γ ) / 2
 x3 = − ( β + γ ) / 2 − p / 3 − 31/2 ( β − γ ) / 2 ; )
( d < 0 ) (  
 δ = ( −d )1/2
 ν = ( b2 − d )1/2
 
( | b | > δ ) (  
 sin(ψ) = δ / ν
 φ = ψ ( b ≥ 0 ) − π ( ( b ≥ 0 ) − 1 ) / 2 ; )
 
( | b | ≤ δ ) (  
 cos(φ) = b / ν ; )
 β = ν1/3 cos(φ/3)
 γ = ν1/3 sin(φ/3)
 x1 = 2 β − p / 3
 x2 = 31/2 γ − β − p / 3
 x3 = − 31/2 γ − β − p / 3 ; )

Se b vale 0 il rapporto b / | b |   vale 1 o −1 ma comunque mai 0, e una disuguaglianza vale 1 se vera e −1 se falsa.

Notare il fatto che la soluzione nel caso di tre radici reali viene ottenuta direttamente, per passaggi matematici, da quella per una sola radice reale.

Viene cioè fatto uso della ben nota relazione tra numeri complessi ossia:

a + b i = ( a2 + b2 )1/2 (cos φ + i sin φ )
cos φ = a / ( a2 + b2 )1/2
sin φ = b / ( a2 + b2 )1/2

per cui la radice cubica di un numero complesso è data da:

( a + b i )1/3= ( a2 + b2 )1/6 (cos φ/3 + i sin φ/3 )

In questo caso nasce il problema numerico che per angoli piccoli l'uso dell' arcocoseno ( Math.acos in Javascript ) è numericamente meno preciso dell'uso dell'arcoseno ( Math.asin in Javascript ) e viceversa, in prossimità di π/2 risulta più preciso l'arcocoseno. Viene dunque effettuata la scelta del metodo migliore il che complica apparentemente ma non sostanzialmente l'algoritmo.

Quando il termine noto è quasi nullo anche una delle radici è prossima a zero ma essendo normalmente ottenuta come differenza di numeri grossi il suo valore potrebbe essere affetto da un errore molto forte in termini di frazione del valore della radice. Un algoritmo solido ma naturalmente più complicato, che tiene conto di questo aspetto è il seguente: 

.
x3 + p x2 + q x + r = 0
p=; q=; r=;

...

Se si ha una sola radice reale si controlla se il suo valore al quadrato risulta inferiore al prodotto delle due radici complesse coniugate e in caso affermativo si calcola la radice piccola come:

x1 = − r / ( x22 + x32 )

dove x2 rappresenta la parte reale e x2 la parte immaginaria di una radice complessa.

Se si hanno viceversa tre radici reali, tranne che nel caso in cui tutte e tre siano nulle, la radice di valore assoluto piu piccola ( supponendo sia x1 ) si ottiene sfruttando il termine noto in modo analogo al caso precedente.

x1 = − r / ( x2 x3 )
dove x2 e x3 rappresentano radici più grandi di x1 in valore assoluto.

Equazione di quarto grado a coefficienti reali

Utilizzo qui la formulazione dell'Abramowitz che mi sembra molto più stabile di altre trovate in giro.

Con questo algoritmo si evita di avere grandezze che tendono all'infinito per particolari valori dei coefficienti dati e dunque i risultati dovrebbero essere affetti da ragionevoli errori numerici. 

.

x4 + a x3 + b x2 + c x + d = 0
a=; b=; c=; d=;

...

Con la cubica trigonometrica...

...

L'algoritmo deriva dai seguenti passi. Supponiamo di voler esprimere il polinomio di quarto grado come prodotto di due trinomi di secondo grado scritti come segue:

x4 + a x3 + b x2 + c x + d = ( x2 + ( p + q ) x + r + s ) ( x2 + ( p − q ) x + r − s )

La scelta dei parametri, un po' ispirata dal cielo, si rivela molto oculata poichè uguagliando i coefficienti dei monomi di ugual grado si trova che deve essere:

a = 2 p
b = 2 r + p2 − q2
c = 2 p r − 2 q s
d = r2 − s2

Ovvero evidenziando meglio i parametri che dobbiamo determinare:

p = a / 2
q2 = 2 r + p2 − b
s2 = r2 − d
q s = p r − c / 2

Il parametro che richiede i calcoli più impegnativi è r che, come si vede, deriva dalla soluzione di una equazione di terzo grado:

( 2 r + p2 − b ) ( r2 − d ) = ( p r − c / 2 )2

ovvero

8 r3 − 4 b r2 + 4 ( c p − 2 d ) r + 4 d ( b − p2 ) − c2 = 0

Per mostrare che la procedura seguita è equivalente (ma numericamente migliore) ad altre, cerchiamo di ottenere l'equazione della risolvente scritta come la si trova in giro. Posto u = r / 2 ed inoltre p = a / 2

u3 − b u2 + ( a c − 4 d ) u + d ( 4 b − a2 ) − c2 = 0

Trovato u e quindi r si nota che il prodotto p s deve dare ( p r − c / 2) e quindi, anche se, noto r, sussiste la libertà di calcolare p ed s indipendentemente, non è lecito scegliere arbitrariamente il segno di entrambe le radici quadrate.

Se ( a r − c ) è una quantità molto piccola, nulla o quasi, è possibile che sia molto prossimo a zero o s o r. Ma la radice quadrata di una quantità molto piccola è affetta da un errore molto più grosso della quantità originale e dunque se, per esempio, s2 ottenuto per differenze di termini grandi vale 1.e-14 mentre il suo vero valore avrebbe dovuto essere 0, si ha che s vale 1.e-7 mentre il vero valore avrebbe dovuto essere 0. In altre parole l'errore è amplificato oltre 10 milioni di volte !

Il metodo da seguire è dunque quello di NON FARE due radici quadrate per dedurre s e q in modo indipendente ma di fare la radice quadrata della quantità più grossa tra le due e dedurre l'altra dividendo ( a r − c ) / 2 per la quantità grossa trovata. In questo modo verranno usati automaticamente anche i segni corretti che si ignorano quando si effettuano entrambe le radici quadrate.

Stesso discorso riguarda la possibilità che s e p non siano coppie di numeri reali ma coppie di immaginari puri. Questa eventualità non può essere esclusa a priori ma la procedura di effettuare una sola radice quadrata e ottenere l'altra tramite divisione permette di trattare correttamente anche questa eventualità.

La possibilità di avere quantità immaginarie nasce quando la risolvente ammette tre radici reali e se ne scelga una tale per cui il quadrato risulti inferiore a 4 d. Ovviamente la strategia da seguire è quella di scegliere sempre la soluzione della risolvente che sia maggiore in valore assoluto ma... la possibilità di trovare soluzioni di equazioni di secondo grado con coefficienti complessi non rende obbligatoria questa scelta ed inoltre evita la necessità di sfruttare la dimostrazione che, in ogni caso, esista una soluzione della risolvente che evita di risolvere equazioni a coefficienti complessi. Ovvio che tale soluzione esiste perché se così non fosse non sarebbe possibile splittare un polinomio di quarto grado a coefficienti reali nel prodotto di due polinomi di secondo grado a coefficienti reali. Ma questo è sempre sicuramente possibile e dunque....

Altra considerazione va fatta quando in termine noto ( ossia d ) è molto piccolo. In questo caso sarebbe meglio non utilizzare la quaterna di radici fornita da questo algoritmo ossia i coefficienti di entrambi i trinomi di secondo grado in cui si fattorizza il polinomio di quarto, ma solo la coppia di radici che ha il modulo più grande.

L' algoritmo dunque ancora più solido di quello prima illustrato effettua automaticamente questa procedura: il termine noto più piccolo in valore assoluto, tra i due trinomi, viene ottenuto come rapporto e non come differenza. Viene inoltre scelta la radice della risolvente di terzo grado che sia la più grande possibile se ce ne sono tre reali. 

.

x4 + a x3 + b x2 + c x + d = 0

a=; b=; c=; d=;

...

Non usando l'algoritmo di Tartaglia...

...

 

Evidenzio qui un grosso problema solo parzialmente superabile: Quando il polinonio è vicino ad una quarta potenza di un polinomio di primo grado, gli algoritmi forniscono risultati PESSIMI ... o CATTIVI !

Il guaio avviene, per esempio con questi dati:

a=12.4 , b=57.66, c=119.164, d=92.3521
per cui una radice vale -3.1

Viceversa cambiando di pochissimo il valore di d la radice cambia drasticamente. Ponendo infatti:

a=12.4 , b=57.66, c=119.164, d=92.352
la radice vale -3.2 per cui cambiando di un milionesimo circa d la radice cambia di un decimo !

Dalle prove fatte mi sembra che ricorrendo alla soluzione dell'equazione cubica con funzioni trigonometriche i risultati siano migliori... anche se resta il fatto che la precisione è molto degradata...

In Mathematichese

Riporto qui come ho programmato queste funzioni nel linguaggio di Mathematica. Ovviamente sarebbe più semplice usare direttamente le funzioni di Mathematica ma... ho incontrato difficoltà pratiche perché Mathematica ( attualmente, versione 9 ) fornisce soluzioni complesse anche quando le radici sono reali ( come per esempio gli autovalori di una matrice SIMMETRICA di ordine tre o quattro ).

Secondo grado

    solve2[r_] := Block[{a = r[[2]], b = r[[1]],de,p,q,s,t,x},
   If[0 >= Abs[b], Return[{0, -a}]];
   de = a^2 - 4 b ;
   p = Re[de];
   q = Im[de];
   If[p > 0, s = Sqrt[(Sqrt[p^2 + q^2] + p)/2]; t = q/(2 s),
    t = Sqrt[(Sqrt[p^2 + q^2] - p)/2];
    If[t == 0, Return[{-a/2, -a/2}]];
    s = q/(2 t)];
   If[s Re[a] + t Im[a] > 0, x = -(a + s + I  t)/2,
    x = -(a - s - I t)/2];
   Return[{x, b /x}]
   ];
poldue[p_, x_] := p[[1]] + x (p[[2]] + x);

Terzo grado trigonometrico

    solve3[poli_] := 
  Block[{a,p,q,A,B,C,D,Ca,s = poli[[3]]/3,sina,cosa,expa, 
    segno, sqrt3 = Sqrt[3]}, 
   q = poli[[1]] - s (poli[[2]]  - s (poli[[3]]  - s)); 
   p = poli[[2]] - s (2  poli[[3]]  - 3 s);
   If[p == 0, a = Sign[q] Abs[ q]^(1/3) - s; Return[{a, a, a}]];
   If[p > 0, B = Sqrt[4 p/3]; D = - q Sqrt[27 /(4 p^3)];
    a = (Log[D + Sqrt[D^2 + 1]])/3;
    expa = Exp[a];
    sina = B (expa - 1/expa)/2;
    cosa = B (expa + 1/expa)/2;
    Return[{sina - s, (-sina +  I sqrt3 cosa)/2 - 
       s, (-sina -  I sqrt3 cosa)/2 - s}]];
   A = Sqrt[-4 p/3];
   C = q Sqrt[-27/(4 p^3)]; Ca = Abs[C];
   If[Ca > 1, segno = C/Ca;
    a = (Log[Ca + Sqrt[Ca^2 - 1]])/3;
    expa = Exp[a];
    sina = segno A (expa - 1/expa)/2;
    cosa = segno  A (expa + 1/expa)/2;
    Return[{- cosa - s,
      (cosa - I sqrt3 sina )/2 - s,
      (cosa + I sqrt3 sina )/2 - s} ]];
   a = ArcSin[C]/3;
   sina = A Sin[a];
   cosa = A Cos[a];
   Return[ {sina - s, (- sina + sqrt3  cosa)/2 - 
      s, (- sina - sqrt3  cosa)/2 - s}]];

Quarto grado trigonometrico

    solve4[p_] := 
 Block[{a = p[[4]],b = p[[3]],c = p[[2]],d = p[[1]],u,s,t,v,w,
    z, e, f, g, h, m, n},
  u = solve3[{4 d b - c^2 - d a^2, c a - 4 d, -b, 1 }];
  z = u[[1]];
  s = a^2/4 + z - b;
  t = z^2/4 - d;
  v = Abs[s]; w = Abs[t];
  If[w > v, w = Sqrt[w]; v = (a z - 2 c)/(4 w),
   If[v > w, t = s; v = Sqrt[v]; w = (a z - 2 c)/(4 v)]];
  If[t >= 0 , e = z/2 - w; f = z/2 + w;
   g = a/2 - v; h = a/2 + v, 
   e = z/2 + I w;
   f = z/2 - I w;
   g = a/2 - I v;
   h = a/2 + I v];
  m = solve2[{e, g}];
  n = solve2[{f, h}];
  Return[{n[[1]], n[[2]], m[[1]], m[[2]]}]
  ]
Ovviamente tutto è ...PERFEZIONABILE ma almeno come esempio di programmazione ...

Ho messo in rete anche un documento CDF ( http://www.elegio.it/cdf/ ) dove sperimento questi algoritmi:

http://www.elegio.it/calcolatrice/quarto-grado-1.zip

Giampaolo Bottoni
gpbottoni@gmail.com
http://www.alumni.polimi.it/it/Wall
( ing. nucleare 1972 )