1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
|
#include "mn_simpson.h"
#include <stdio.h>
//************************************************************
// FUNCION PARA CALCULAR UNA INTEGRAL POR EL METODO DE SIMPSON
//************************************************************
real mn_simpson(
real (*f)(const real x) /* funcion que se integra */,
const real a /* extremo izquierdo del intervalo */,
const real b /* extremo derecho del invervalo*/,
const int N /* numero de intervalos para calcular la integral */)
{
real integral=0; // variable donde se acumula el valor de la integral
real f_k=(*f)(a); // evaluacion de la funcion en a
real h=(b-a)/(6.*N); // factor multiplicativo
for(int i=0;i<N;i++) // bucle para acumular el valor de la integral
{
real xk=a+i*(b-a)/N; // extremo izquierdo subintervalo
real xk1=a+(i+1)*(b-a)/N; // extremo derecho subintervalo
real f_k_1=(*f)(xk1); // evaluación de la función en el extremo derecho
integral+=h*(f_k+f_k_1+4*(*f)((xk+xk1)/2.)); // acumulacion valor integral
f_k=f_k_1; // actualizacion valor integral extremo izquierdo
}
return integral; // se devuelve el valor de la integral
}
//*************************************************************************
// FUNCION PARA CALCULAR UNA INTEGRAL POR ITERACIONES DEL METODO DE SIMPSON
//*************************************************************************
real mn_simpson(
real (*f)(const real x) /* funcion que se integra */,
const real a /* extremo izquierdo del intervalo */,
const real b /* extremo derecho del invervalo*/,
int &Niter /* variable de salida con el numero de iteraciones realizadas */,
const real TOL /* tolerancia para controlar la convergencia */)
{
int N=1;
real integral =mn_simpson((*f),a,b,N);
real error=TOL+1;
Niter=0;
while(error>TOL ){
if(N>1e7) break; // fijamos "a priori" el numero máximo de intervalo a 1e7
Niter++;
N*=2;
real integral2=mn_simpson((*f),a,b,N);
error=mn_distancia(integral,integral2);
integral=integral2;
//printf("integral=%e error=%e N=%d \n",(double) integral,(double) error,N);
}
return(integral);
}
|
