# 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_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_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