Last Update: Feb. 9, 2006
#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;
}
$ ./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
#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;
}
$ ./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