newton_s_method.c

来自「A program for solving equations using Ne」· C语言 代码 · 共 103 行

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#define TOL 10.0e-6

void deri(double coef[], int deg)
{
    int i;
    for (i=1; i<=deg; i++)
    {
        coef[i-1] = i * coef[i];
    }
    coef[i-1] = 0.0;
}

double poly(double coef[], int deg, double x)
{
    int k = 0;
    double y = coef[deg];
    for(k=deg; k>0; k--)
    {
        y = y*x + coef[k-1];
    }
    return y;
}

double Newton(double f[], int deg, double p0)
{
    int j;
    double fx = 0.0;
    double dfx = 0.0;
    double df[deg+1];
    for(j=deg; j>=0; j--)
    {
        df[j] = f[j];
    }   
    deri(&df, deg);
    printf("p0 = %.8lf\n", p0);
    j=1;
    while( poly(f, deg, p0) > TOL || poly(f, deg, p0) < TOL*(-1.0))
    {
        fx = poly(f, deg, p0);
        dfx = poly(df, deg-1, p0);
       
        if(dfx != 0.0)
            p0 = p0 - fx/dfx;
        else
        {
            printf("denormitor error!\n");
            break;
        }
        printf("p%d = %.8lf\n", j++, p0);
    }
    printf("f(p) = %f\n", poly(f, deg, p0));
    return p0;
}

double Secant(double f[], int deg, double p0, double p1)
{
    int j;
    double fx0, fx1;
    double tmp = p1;
    double df[deg+1];
    for(j=deg; j>=0; j--)
    {
        df[j] = f[j];
    }   
    deri(&df, deg);
    printf("p0 = %.8lf\n", p0);
    printf("p1 = %.8lf\n", p1);
    j=2;
    while(poly(f, deg, p1) > TOL || poly(f, deg, p1) < TOL*(-1.0))
    {
        fx0 = poly(f, deg, p0);
        fx1 = poly(f, deg, p1);
        if(fx1-fx0 != 0.0)
        {
            tmp = p1;
            p1 = p1 - fx1*(p1-p0)/(fx1-fx0);
            p0 = tmp;
        }
        else
        {
            printf("denormitor error!\n");
            break;
        }
        printf("p%d = %.8lf\n", j++, p1);
    }
    printf("f(p) = %.8lf\n", poly(f, deg, p1));
    return p1;
}

int main()
{
    double test[] = {-9.0, -221.0, +9.0, +18.0, 230.0};
//    double test_1[] = {1.0, 0.0, 0.0, -0.5};
//    printf("%.8lf\n", Newton(test, 4,  -0.5));
//    printf("%.8lf\n", Newton(test, 4,   0.5));
//    printf("%.8lf\n", Secant(test, 4, -1.0, 0.0));
//    printf("%.8lf\n", Secant(test, 4,  0.0, 1.0));
}