[C, Algoritmi] Metodo di Newton e bisezione

[lucia88]

Ciao a tutti,
ho il seguente esercizio
1) Si consideri f(x)=e^(-1)+ x^2 -2:
-b: scrivere un algoritmo che stimi la soluzione con almeno due metodi

Ho provato a scrivere il seguente programma che calcola una stima delle radici con il metodo di newton e il metodo di bisezione, però quando inserisco gli intervalli I[-1,0] e B[1,2] (che sono gli intervalli che contengono una radice) vengono dei risultati strani :
Ho sbagliato l'algoritmo di newton o bisezione secondo voi?

#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<iostream>
using namespace std;
 
const double e=2.7182;
 // funzione
double f(double x)
{
                return (pow(e,-x)+pow(x,2)-2);
}
 
 
//derivata prima
double f1(double x)
{
                return (pow(-e,-x)+2*x);
}
 
//derivata secondA
double f2(double x)
{
                return(pow(e,-x)+2);
}
 
//funzione che legge la precisione,cioè errore massimo, e il numero max d'iterazioni
double precisione(int& n)
{
                double eps;
                do
                {
                               cout<<endl<<"Inserire il numero massimo d'iterazioni:";
                               cin>>n;
                }while(n<=0);
               
                do
                {
                    cout<<endl<<"Inserire la precisione(deve essere maggiore di zero):";
                    cin>>eps;       
                }while(eps<=0);
               
                return eps;
}
 
//funzione che legge gli estremi dell'intervallo in cui cercare lo zero e restituisce il p.innesco
//per metodo di newton
double range(double& a, double& b)
{
                cout<<endl<<"Inserire un intervallo [a,b] all'interno del quale si trova lo zero,"<<endl;
                cout<<"e nel quale l'algoritmo per la ricerca degli zeri sia convergente:";
                cout<<endl;
               
                do
                {
                               cout<<"a=";
                               cin>>a;
                              
                               cout<<"b=";
                               cin>>b;
                              
                               if(f(a)*f(b)>0)
                            cout<<"Attenzione: f(a) ha lo stesso segno di f(b).Riprovare...";
                }while(f(a)*f(b)>0);
               
                if(f(a)*f2(a)>0) return b;
                else return a;
}
 
void stampa(int n,double x,double err)
{
                cout<<endl<<"Radice calcolata:"<<endl<<x;
                cout<<endl<<"Numero d'iterazioni:"<<endl<<n;
                cout<<endl<<"Errore stimato(Tolleranza sul test d'arresto):"<<endl<<err;
                return;
}
 
double bisezione(double a, double b,int nmax,double err,int& cont)
{
                cont=0;
                do
                {
                               if(fabs(f(a))/f(a)==fabs(f((a+b)/2))/f((a+b)/2)) a=(a+b)/2;
                               else b=(a+b)/2;
                               cont++;
                               //criterio d'arresto: controllo del residuo
                }while(cont<nmax&&fabs(f((a+b)/2))>err);
               
                return (a+b)/2;
}
 
int newton(double a,double err,int nmax,double& x)
{
                int cont=0;
                double tmp;
               
                do
                {
                               if(f1(a)==0)
                               {
                                               cout<<endl<<"errore:il programma verra' terminato...";
                                               exit(1);
                               }
                              
                               x=a+f(a)/f1(a);
                               cont++;
                               tmp=a;
                               a=x;
                               x=tmp;
                //criterio d'arresto: controllo dell'incremento
                }while(cont<nmax&&fabs(a-x)>err);
               
                return cont;
}
 
int main()
{
                //nmax= numero massimo d'iterazioni
                //it_N=iterazioni del metodo di newton
                //it_B=iterazioni del metodo di bisezione
                int nmax,it_N,it_B;
               
                //toll=errore massimo
                //a,b=estremi intervallo
                //x0=punto di innesco del metodo di newton
                //risultato del metodo di newton
                double a,b,toll,x0,x;
               
                //leggo da testiera precisione e numero massimo d'iterazioni
               
                toll=precisione(nmax);
               
                //leggo l'intervallo e il punto di innesco di newton:
               
                x0=range(a,b);
               
                //calcolo la radice con newton,memorizzando anche il numero di iterazioni svolte
               
                it_N=newton(x0,toll,nmax,x);
               
                //calcolo la radice con il metodo di bisezione, memorizzo anche il numero di iterazioni
               
                b=bisezione(a,b,nmax,toll,it_B);
               
                //stampo su schermo i risultati ottenuti con i due metodi:
               
                cout<<endl<<"METODO DI NEWTON:";
                cout<<endl;
                stampa(it_N,x,toll);
               
                cout<<endl<<"METODO DI BISEZIONE:";
                cout<<endl;
                stampa(it_B,b,toll);
               
                return 0;
}

[claudio862]

Innanzitutto hai sbagliato la funzione che calcola la derivata prima:

double f1(double x)
{
    // Guarda dove è il segno -
    //return (pow(-e, -x) + 2*x);
    return (-pow(e, -x) + 2*x);
}

Poi hai sbagliato la funzione che calcola la radice con il metodo di bisezione. Probabilmente ti aspetti che fabs(x)/x dia +1.0 oppure -1.0 a seconda se x sia positivo o meno. Purtroppo le divisioni con i floating point danno errori di approssimazione, quindi la tua condizione nell'if:

if(fabs(f(a))/f(a)==fabs(f((a+b)/2))/f((a+b)/2))
...

sarà sempre falsa, anche se hanno lo stesso segno, perché i due risultati potrebbero essere ad esempio 1.0000000003 e 1.000000000005. Infatti controllare se due variabili floating point sono uguali è quasi sempre un errore.

Il modo classico per vedere se due numeri hanno lo stesso segno è di moltiplicarli e vedere se il risultato è positivo.

double bisezione(double a, double b, int nmax, double err, int & cont)
{
    double midpoint, midpointValue;
    cont = 0;
    do
    {
        midpoint = (a + b) / 2;
        const double leftValue = f(a);
        midpointValue = f(midpoint);

        if (leftValue * midpointValue > 0) {
            // f((a+b)/2) ha lo stesso segno di f(a), quindi
            // la soluzione è tra (a+b)/2 e b.
            a = midpoint;
        } else {
            // f((a+b)/2) ha segno diverso da f(a), quindi
            // esiste una soluzione tra a e (a+b)/2.
            b = midpoint;
        }
        cont++;

        //criterio d'arresto: controllo del residuo
    } while (cont < nmax && fabs(midpointValue) > err);

    return midpoint;
}

Infine, hai sbagliato la funzione che calcola la radice con il metodo di Newton Il passo di aggiornamento è più o meno giusto, la condizione di arresto no, e il ruolo di "a" mi sembra abbastanza contorto (e comunque sbagliato).

int newton(double a, double err, int nmax, double & x)
{
    int cont = 0;

    x = a; // a è la soluzione iniziale x0 del metodo.

    do
    {
        const double derivata = f1(x);

        if (derivata == 0) // Stai eguagliando due double, non saranno quasi mai uguali, come detto prima.
        {
            cout << "errore:il programma verra' terminato..." << '\n';
            exit(1);
        }

        x = x - f(x) / derivata;
        cont++;

    //criterio d'arresto: controllo dell'incremento
    } while (cont < nmax && fabs(f(x)) > err);

    return cont;
}

Poi qualche commento sul codice:

- Non indentare così tanto. Stai sui 4 o 8 spazi per ogni livello.

- Usa di più gli spazi nelle espressioni, non scrivere così

if(fabs(f(a))/f(a)==fabs(f((a+b)/2))/f((a+b)/2)) ...

ma così

if (fabs(f(a)) / f(a) == fabs(f((a + b) / 2)) / f((a + b) / 2)) ...

- Non scrivere espressioni troppo complicate, usa variabili di appoggio:

const double segnoPuntoSinistro = fabs(f(a)) / f(a);
const double segnoPuntoMediano = fabs(f((a + b) / 2)) / f((a + b) / 2);
if (segnoPuntoSinistro == segnoPuntoMediano) ...

(fermo restando che questo if è sbagliato, come detto sopra).

- Non andare a capo prima di scrivere qualcosa, fallo dopo. Quindi, non scrivere così

void stampa(int n,double x,double err)
{
    cout << endl << "Radice calcolata:" << endl << x;
    cout << endl << "Numero d'iterazioni:" << endl << n;
    cout << endl << "Errore stimato:" << endl << err;
}

ma così

void stampa(int n, double x, double err)
{
    cout << "Radice calcolata:"    << endl << x    << endl;
    cout << "Numero d'iterazioni:" << endl << n    << endl;
    cout << "Errore stimato:"      << endl << err  << endl;
    cout << "Errore vero:"         << endl << f(x) << endl;
}

- Cerca di mantenere un'interfaccia comune dove possibile. Hai due funzioni che calcolano la radice di una funzione. Il metodo di Newton restituisce il numero di iterazioni e scrive la radice in una variabile passata come riferimento, mentre il metodo di bisezione fa il contrario.

- Non dichiarare tutte le variabili all'inizio di una funzione. Fallo quando ti servono (e se non variano dichiarale come costanti).

[vict85]

Il metodo di Claudio per il segno ha il lato negativo che potrebbe essere affetto da overflow o underflow. Potrebbe quindi essere comunque opportuno fare qualcosa del tipo \(a\lvert a\rvert^{-1}b\). D'altra parte dipende dal criterio di arresto usato: se err è sufficientemente grande non c'è pericolo di underflow. Il pericolo di overflow viene invece ignorato per il fatto di trovarci vicino a zero. Vorrei comunque far notare che la convergenza del metodo di newton depende dall'intervallo considerato. È quindi alle volte utile fare qualche passo con qualche altro metodo prima di usarlo. Ma sono questioni di analisi numerica.[strike][/strike]

Detto questo esiste la funzione exp ed è meglio usare quella che pow con un valore di $e$ molto arrotondato. Tieni inoltre conto che in pow se la base è negativa e l'esponente non intero il risultato è NaN.

[lucia88]

Grazie mille per le correzioni e i suggerimenti
L'ho corretto e adesso il programma funziona meglio...

Ho poi aggiunto al programma delle istruzioni per visualizzare il grafico della funzione con GNUplot, lanciandolo con l'istruzione:

system("start C:\\Users\\lucia\\Deskop\\wgnuplot.exe comandi.txt");

(dove il file comandi.txt contiene le istruzioni per gnuplot)

però quando eseguo il programma il grafico di gnuplot appare e scompare immediatamente,
come faccio per ''fissarlo sullo schermo''?