diff options
Diffstat (limited to 'doc/libinum-examples/src')
| -rw-r--r-- | doc/libinum-examples/src/integral_gauss-legendre.c | 92 | ||||
| -rw-r--r-- | doc/libinum-examples/src/root_newton.c | 58 | 
2 files changed, 150 insertions, 0 deletions
diff --git a/doc/libinum-examples/src/integral_gauss-legendre.c b/doc/libinum-examples/src/integral_gauss-legendre.c new file mode 100644 index 0000000..f65465f --- /dev/null +++ b/doc/libinum-examples/src/integral_gauss-legendre.c @@ -0,0 +1,92 @@ +#include <stdio.h> +#include <stdarg.h> +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include <mpfr.h> +#include <libinum.h> + +int func(mpfr_t* out, mpfr_t in, void* extra_args); + +int main(int argc, const char* argv[]){ +  int ret; +  mpfr_t tolerance, val, lower, upper; +  array_mpfr abcissa, weights; + +  // precision of MPFR floats +  mpfr_set_default_prec(53); + +  // tolerance for computing weights +  mpfr_init(tolerance); +  mpfr_set_d(tolerance, 1.e-11, MPFR_RNDN); + +  // compute weights +  ret=gauss_legendre_weights_mpfr(10, tolerance, 10000, &abcissa, &weights); +  // return codes +  if(ret==LIBINUM_ERROR_MAXITER){ +    printf("error: maximum number of iterations reached when computing the integration abcissa\n"); +    mpfr_clear(tolerance); +    return(ret); +  } +  else if(ret==LIBINUM_ERROR_NAN){ +    printf("error: infinity encountered when computing the integration abcissa\n"); +    mpfr_clear(tolerance); +    return(ret); +  } +  else{ +    printf("abcissa:\n"); +    array_mpfr_print(abcissa); +    printf("\nweights:\n"); +    array_mpfr_print(weights); +    printf("\n"); +  } + +  mpfr_clear(tolerance); + +  mpfr_inits(val, lower, upper, NULL); + +  mpfr_set_ui(lower, 0, MPFR_RNDN); +  mpfr_set_ui(upper, 2, MPFR_RNDN); + +  ret=integrate_gauss_mpfr(&val, &func, lower, upper, abcissa, weights, NULL); +  // return codes +  if(ret==LIBINUM_ERROR_SIZE_MISMATCH){ +    printf("error: the mpfr_arrays 'abcissa' and 'weights' must have the same length\n"); +    mpfr_clears(val, upper, lower, NULL); +    array_mpfr_free(abcissa); +    array_mpfr_free(weights); +    return(ret); +  } +  else if(ret==LIBINUM_ERROR_NAN){ +    printf("error: the integrand is singular\n"); +    array_mpfr_free(abcissa); +    array_mpfr_free(weights); +    mpfr_clears(val, upper, lower, NULL); +    return(ret); +  } +  else{ +    printf("\\int_0^2 x/sin(x) dx = "); +    fprint_mpfr(stdout, val); +    printf("\n"); +  } + +  mpfr_clears(val, upper, lower, NULL); +  array_mpfr_free(abcissa); +  array_mpfr_free(weights); + +  return(0); +} + +/* +// e^x +int func(mpfr_t* out, mpfr_t in, void* extra_args){ +  mpfr_exp(*out, in, MPFR_RNDN); +  return(0); +} +*/ + +// x/sin(x) +int func(mpfr_t* out, mpfr_t in, void* extra_args){ +  mpfr_sin(*out, in, MPFR_RNDN); +  mpfr_div(*out, in, *out, MPFR_RNDN); +  return(0); +} diff --git a/doc/libinum-examples/src/root_newton.c b/doc/libinum-examples/src/root_newton.c new file mode 100644 index 0000000..f307c58 --- /dev/null +++ b/doc/libinum-examples/src/root_newton.c @@ -0,0 +1,58 @@ +#include <stdio.h> +#include <stdarg.h> +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include <mpfr.h> +#include <libinum.h> + +int func(mpfr_t* out, mpfr_t in, void* extra_args); +int dfunc(mpfr_t* out, mpfr_t in, void* extra_args); + +int main(int argc, const char* argv[]){ +  int ret; +  mpfr_t init, prec, x; +  unsigned int maxiter; +  int extra_arg; + +  mpfr_inits(init, prec, x, NULL); + +  // start at 2 +  mpfr_set_ui(init, 2, MPFR_RNDN); +  // precision +  mpfr_set_d(prec, 1.e-30, MPFR_RNDN); +  // maximum number of iterations before error +  maxiter=10000; +  // extra argument +  extra_arg=4; + +  // compute root +  ret=root_newton_mpfr(&x, &func, &dfunc, init, prec, maxiter, &extra_arg); + +  // return codes +  if(ret==LIBINUM_ERROR_MAXITER){ +    printf("error: maximum number of iterations reached\n"); +  } +  else if(ret==LIBINUM_ERROR_NAN){ +    printf("error: infinity encountered\n"); +  } +  else{ +    mpfr_printf("% 14.7Re\n", x); +  } + +  mpfr_clears(init, prec, x, NULL); + +  return(0); +} + +// x^2-1 +int func(mpfr_t* out, mpfr_t in, void* extra_args){ +  mpfr_pow_ui(*out, in, 2, MPFR_RNDN); +  mpfr_add_d(*out, *out, -1./ *((int*)extra_args), MPFR_RNDN); +  return(0); +} + +// 2*x +int dfunc(mpfr_t* out, mpfr_t in, void* extra_args){ +  mpfr_mul_ui(*out, in, 2, MPFR_RNDN); +  return(0); +}  | 
