Solving a nonlinear equation with Newton method

Last Update: Feb. 9, 2006


MPFR

Sample code

newton_mpf.c

#include <stdio.h>
#include <math.h>

/* Header file of BNCpack */
#include "bnc.h"

/* solve f(x) = 0 */

/* f(x) */
void mpffunc(MPFVector vret, MPFVector x)
{
    mpf_t tmp_x[3], tmp[3];

    /* Initialize */
    mpf_init(tmp_x[0]);
    mpf_init(tmp_x[1]);
    mpf_init(tmp_x[2]);
    mpf_init(tmp[0]);
    mpf_init(tmp[1]);
    mpf_init(tmp[2]);

    mpf_set(tmp_x[0], get_mpfvector_i(x, 0));
    mpf_set(tmp_x[1], get_mpfvector_i(x, 1));
    mpf_set(tmp_x[2], get_mpfvector_i(x, 2));

    /* f_1(x) = x1 + x2 + x3 - 6        */
    /* f_2(x) = x1 * x2 * x3 - 6        */
    /* f_3(x) = x1^2 + x2^2 + x3^2 - 14 */

    mpf_set(tmp[0], tmp_x[0]);
    mpf_add(tmp[0], tmp[0], tmp_x[1]);
    mpf_add(tmp[0], tmp[0], tmp_x[2]);
    mpf_sub_ui(tmp[0], tmp[0], 6UL);
    set_mpfvector_i(vret, 0, tmp[0]);

    mpf_set(tmp[1], tmp_x[0]);
    mpf_mul(tmp[1], tmp[1], tmp_x[1]);
    mpf_mul(tmp[1], tmp[1], tmp_x[2]);
    mpf_sub_ui(tmp[1], tmp[1], 6UL);
    set_mpfvector_i(vret, 1, tmp[1]);

    mpf_mul(tmp[0], tmp_x[0], tmp_x[0]);
    mpf_mul(tmp[1], tmp_x[1], tmp_x[1]);
    mpf_mul(tmp[2], tmp_x[2], tmp_x[2]);
    mpf_add(tmp[2], tmp[2], tmp[0]);
    mpf_add(tmp[2], tmp[2], tmp[1]);
    mpf_sub_ui(tmp[2], tmp[2], 14UL);
    set_mpfvector_i(vret, 2, tmp[2]);

    /* Clean */
    mpf_clear(tmp_x[0]);
    mpf_clear(tmp_x[1]);
    mpf_clear(tmp_x[2]);
    mpf_clear(tmp[0]);
    mpf_clear(tmp[1]);
    mpf_clear(tmp[2]);
}

/* jacobian matrix of f(x) */
void jac_mpffunc(MPFMatrix mret, MPFVector x)
{
    mpf_t tmp_x[3], tmp;

    /* Initialize */
    mpf_init(tmp_x[0]);
    mpf_init(tmp_x[1]);
    mpf_init(tmp_x[2]);
    mpf_init(tmp);

    mpf_set(tmp_x[0], get_mpfvector_i(x, 0));
    mpf_set(tmp_x[1], get_mpfvector_i(x, 1));
    mpf_set(tmp_x[2], get_mpfvector_i(x, 2));

    /* [       1       1       1 ] */
    /* [ x2 * x3 x1 * x3 x1 * x2 ] */
    /* [  2 * x1  2 * x2  2 * x3 ] */

    set_mpfmatrix_ij_ui(mret, 0, 0, 1UL);
    set_mpfmatrix_ij_ui(mret, 0, 1, 1UL);
    set_mpfmatrix_ij_ui(mret, 0, 2, 1UL);

    mpf_mul(tmp, tmp_x[1], tmp_x[2]);
    set_mpfmatrix_ij(mret, 1, 0, tmp);

    mpf_mul(tmp, tmp_x[0], tmp_x[2]);
    set_mpfmatrix_ij(mret, 1, 1, tmp);

    mpf_mul(tmp, tmp_x[0], tmp_x[1]);
    set_mpfmatrix_ij(mret, 1, 2, tmp);

    mpf_mul_ui(tmp, tmp_x[0], 2UL);
    set_mpfmatrix_ij(mret, 2, 0, tmp);

    mpf_mul_ui(tmp, tmp_x[1], 2UL);
    set_mpfmatrix_ij(mret, 2, 1, tmp);

    mpf_mul_ui(tmp, tmp_x[2], 2UL);
    set_mpfmatrix_ij(mret, 2, 2, tmp);

    /* Clean */
    mpf_clear(tmp_x[0]);
    mpf_clear(tmp_x[1]);
    mpf_clear(tmp_x[2]);
    mpf_clear(tmp);
}

/* Main function */
int main()
{
    long int itimes;
    mpf_t mpfabs_tol, mpfrel_tol;
    MPFVector mpfx_init, mpfx_ans;

    /* Set default precision in decimal */
    set_bnc_default_prec_decimal(50);

    /* Initialize */
    mpfx_init = init_mpfvector(3);
    mpfx_ans  = init_mpfvector(3);
    mpf_init(mpfabs_tol);
    mpf_init(mpfrel_tol);

    /* Set initial guess */
    set_mpfvector_i_ui(mpfx_init, 0, 6UL);
    set_mpfvector_i_ui(mpfx_init, 1, 5UL);
    set_mpfvector_i_ui(mpfx_init, 2, 4UL);

    /* Print initial guess */
    printf("Initial guess:\n");
    print_mpfvector(mpfx_init);

    /* Set absolute and relative tolerances */
    mpf_set_d(mpfabs_tol, 1.0e-200);
    mpf_set_d(mpfrel_tol, 1.0e-40);

    /* call Newton Method */
    itimes = mpf_newton(
        mpfx_ans,   /* One of approximate zeros */
        mpfx_init,  /* Initial guess */
        mpffunc,    /* f(x) */
        jac_mpffunc,/* Jacobian matrix function of f(x) */
        100,        /* Maximum iterative times */
        mpfabs_tol, /* Absolute tolerance */
        mpfrel_tol  /* Relative tolerance */
    );

    /* Print mpfx_ans */
    printf("Approximate zero(Iterative times %d): \n", itimes);
    print_mpfvector(mpfx_ans);

    /* Clear */
    free_mpfvector(mpfx_init);
    free_mpfvector(mpfx_ans);
    mpf_clear(mpfabs_tol);
    mpf_clear(mpfrel_tol);

    /* Exit */
    return 0;
}

Example execution

$ ./newton_mpf
-------------------------------------------------------------------------------
BNC Default Precision    : 167 bits(50.3 decimal digits)
BNC Default Rounding Mode: Round to Nearest
-------------------------------------------------------------------------------
Initial guess:
    0 6.000000000000000000000000000000000000000000000000000
    1 5.000000000000000000000000000000000000000000000000000
    2 4.000000000000000000000000000000000000000000000000000
Approximate zero(Iterative times 19): 
    0 1.999999999999999999999999999999999999999999999999979
    1 1.000000000000000000000000000000000000000000000000011
    2 3.000000000000000000000000000000000000000000000000000

IEEE754 double

Sample code

newton_d.c

#include <stdio.h>
#include <math.h>

/* Header file of BNCpack */
#include "bnc.h"

/* solve f(x) = 0 */

/* f(x) */
void dfunc(DVector vret, DVector x)
{
    double tmp_x[3];

    tmp_x[0] = get_dvector_i(x, 0);
    tmp_x[1] = get_dvector_i(x, 1);
    tmp_x[2] = get_dvector_i(x, 2);

    /* f_1(x) = x1 + x2 + x3 - 6        */
    /* f_2(x) = x1 * x2 * x3 - 6        */
    /* f_3(x) = x1^2 + x2^2 + x3^2 - 14 */

    set_dvector_i(vret, 0, tmp_x[0] + tmp_x[1] + tmp_x[2] - 6.0);
    set_dvector_i(vret, 1, tmp_x[0] * tmp_x[1] * tmp_x[2] - 6.0);
    set_dvector_i(vret, 2, tmp_x[0] * tmp_x[0] + tmp_x[1] * tmp_x[1] + tmp_x[2] * tmp_x[2] - 14.0);
}

/* jacobian matrix of f(x) */
void jac_dfunc(DMatrix mret, DVector x)
{
    double tmp_x[3];

    tmp_x[0] = get_dvector_i(x, 0);
    tmp_x[1] = get_dvector_i(x, 1);
    tmp_x[2] = get_dvector_i(x, 2);

    /* [       1       1       1 ] */
    /* [ x2 * x3 x1 * x3 x1 * x2 ] */
    /* [  2 * x1  2 * x2  2 * x3 ] */

    set_dmatrix_ij(mret, 0, 0, 1.0);
    set_dmatrix_ij(mret, 0, 1, 1.0);
    set_dmatrix_ij(mret, 0, 2, 1.0);

    set_dmatrix_ij(mret, 1, 0, tmp_x[1] * tmp_x[2]);
    set_dmatrix_ij(mret, 1, 1, tmp_x[0] * tmp_x[2]);
    set_dmatrix_ij(mret, 1, 2, tmp_x[0] * tmp_x[1]);

    set_dmatrix_ij(mret, 2, 0, 2.0 * tmp_x[0]);
    set_dmatrix_ij(mret, 2, 1, 2.0 * tmp_x[1]);
    set_dmatrix_ij(mret, 2, 2, 2.0 * tmp_x[2]);
}

/* Main function */
int main()
{
    long int itimes;
    double dabs_tol, drel_tol;
    DVector dx_init, dx_ans;

    /* Initialize */
    dx_init = init_dvector(3);
    dx_ans  = init_dvector(3);

    /* Set initial guess */
    set_dvector_i(dx_init, 0, 6.0);
    set_dvector_i(dx_init, 1, 5.0);
    set_dvector_i(dx_init, 2, 4.0);

    /* Print initial guess */
    printf("Initial guess:\n");
    print_dvector(dx_init);

    /* Set absolute and relative tolerances */
    dabs_tol = 1.0e-50;
    drel_tol = 1.0e-8;

    /* call Newton Method */
    itimes = dnewton(
        dx_ans,     /* One of approximate zeros */
        dx_init,    /* Initial guess */
        dfunc,      /* f(x) */
        jac_dfunc,  /* Jacobian matrix function of f(x) */
        100,        /* Maximum iterative times */
        dabs_tol,   /* Absolute tolerance */
        drel_tol    /* Relative tolerance */
    );

    /* Print dx_ans */
    printf("Approximate zero(Iterative times %d): \n", itimes);
    print_dvector(dx_ans);

    /* Clear */
    free_dvector(dx_init);
    free_dvector(dx_ans);

    /* Exit */
    return 0;
}

Example execution

$ ./newton_d
Initial guess:
    0   6.00000000000000000e+00
    1   5.00000000000000000e+00
    2   4.00000000000000000e+00
Approximate zero(Iterative times 17): 
    0   1.99999999999999445e+00
    1   1.00000000000000089e+00
    2   3.00000000000000444e+00

<- Go back