viernes, 5 de mayo de 2017

Integración Numérica


Antecedentes
La integral definida de una función, partiendo de un intervalo, es igual al área que se encuentra debajo de la curva hasta el eje x hasta el extremo del intervalo. Los primeros métodos para calcularlo consisten en dividir la selección en figuras geométricas donde es más sencillo calcular su área. Dado que las áreas debajo de cualquier función no tienen líneas rectas y mayoritariamente son tienen distintas tipos de curvas, este tipo de métodos solo calculan aproximaciones.
La suma de Riemann, por ejemplo, consiste en formar rectángulos a lo largo del intervalo dado y calcular el área de cada una de ellas. La altura de estos rectángulos puede ser considerada como el primer o segundo punto evaluado en la función. He ésta, por ejemplo, la diferencia entre los métodos de suma de rectángulos a la izquierda y derecha.
Justificación y Propósito
Podemos pensar en el área debajo de una curva como la suma de muchos pequeños rectángulos. ¿Qué tan pequeños deben de ser estos rectángulos? Depende...


Así es como se tomarían los rectángulos a la izquierda:
Y así es como se tomarían los rectángulos a la derecha:
Podría, entonces, creerse que entre más pequeños sean los rectángulos, más exacto será nuestra aproximación. Sin embargo, esto no siempre es cierto (como se observará en los ejemplos más adelante). Esta intuición sí es verdadera para otro método de suma… en vez de sumar rectángulos, y sobreestimar o subestimar parte de estos rectángulos, ¿por qué no usamos trapecios? Esto reduce nuestro error:
Diagrama de flujo

Suma de Riemann

Riemann.png

Regla del Trapecio
Trapezoidal.png
Código en C++

 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
52
53
54
55
56
57
58
59
60
61
62
63
64
#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;

double f(double x){
  // return 0.2 + (25*x) - (200*pow(x,2)) + (675*pow(x,3)) - (900*pow(x,4)) + (400*pow(x,5));
  return exp(pow(x,2));
}

double areaRectLeft(double a, double b, double div){
  double area = 0.0;
  double rectangleSize = (b-a)/div;
  b = a+rectangleSize;
  for(int i=0; i<div; i++){
    area += ( (f(a))*(b-a) ); // area of each rectangle
    // cout << "a:" << a  << " b:" << b << endl;
    a += rectangleSize;
    b += rectangleSize;
  }
  return area;
}

double areaRectRight(double a, double b, double div){
  double area = 0.0;
  double rectangleSize = (b-a)/div;
  b = a+rectangleSize;
  for(int i=0; i<div; i++){
    area += ( (f(b))*(b-a) ); // area of each rectangle
    // cout << "a:" << a  << " b:" << b << endl;
    a += rectangleSize;
    b += rectangleSize;
  }
  return area;
}

double areaTrap(double a, double b, double div){
  double area = 0.0;
  double trapSize = (b-a)/div;
  b = a+trapSize;
  for(int i=0; i<div; i++){
    area += ( (f(a)+f(b))*(b-a) )/2; // area of each trapezoid
    // cout << "a:" << a  << " b:" << b << endl;
    a += trapSize;
    b += trapSize;
  }
  return area;
}

int main(){
  cout << endl << "RECTANGULAR IZQUIERDA" << endl;
  for(int i=1; i<16; i++){
    cout << i << ": " << "area under curve: " << areaRectLeft(0,1,i) << "     " << 1.46265-areaRectLeft(0,1,i) << endl;
  }
  cout << endl << "RECTANGULAR DERECHA" << endl;
  for(int i=1; i<16; i++){
    cout << i << ": " << "area under curve: " << areaRectRight(0,1,i) << "     " << 1.46265-areaRectRight(0,1,i) << endl;
  }
  cout << endl << "TRAPECIO" << endl;
  for(int i=1; i<16; i++){
    cout << i << ": " << "area under curve: " << areaTrap(0,1,i) << "     " << 1.46265-areaTrap(0,1,i) << endl;
  }
  return 0;
}



 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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

// Method that evaluates the x in a function
double f( double x  ) {

    //return 0.2 + 25 * x - 200 * pow( x, 2 ) + 675 * pow( x, 3 ) - 900 * pow( x, 4 ) + 400 * pow( x, 5 );
    return exp ( x * x );
}

void SimpsonRule( int n, double *x ) {

    cout << "*** Simpson Rule 1/3 ***" << endl;

    double l = 0;
    double oddSum = 0;
    double pairSum = 0;

    for ( int i = 1; i < n - 1; i++ ) {

        if ( i % 2 == 0 )
            pairSum += f( x[i] );
        else
            oddSum += f( x[i] );

    }

    l = ( x[n - 1] - x[0] ) * ( f( x[0] ) + 4 * oddSum + 2 * pairSum + f( x[n - 1] ) ) / ( 3 * (n - 1) );

    cout << "Area under the function = " << l << endl;
}

void SimpsonRule3() {

    double x[4] = {0, 0, 0, 0};
    int n = 4;

    cout << "*** Simpson Rule 3/8 ***" << endl;

    double l = 0;

    l = ( x[n - 1] - x[0] ) * ( f(x[0]) + 3 * f(x[1]) + 3 * f(x[2]) + f(x[3]) ) / 8;

    cout << "Area under the function = " << l << endl;
}

int main() {

    //double x[4] = {0, 0.5, 1};

    //int n = sizeof(x) / sizeof(x[0]);

    double *x;

    double a = 0;
    double b = 1;

    int n = 4;  //<- Value to be changed

    x = new double[n + 1];

    double frag = ( b - a ) / n;
    for ( int i = 1; i < n; i++ ) {

        x[i] = x[i - 1] + frag;
    }

    x[0] = a;
    x[n] = b;

    n++;

    SimpsonRule( n, x );
    return 0;
}


Ejemplo Comparativo
Encontremos, por ejemplo, el área bajo la siguiente desde 0 a 1 de la siguiente curva:




Área real: 1.46265
Veamos a los resultados de diferentes métodos con diferente cantidad de divisiones:
divisiones
rectangular izquierda
error absoluto
rectangular derecha
error absoluto
trapecio
error absoluto
1
1
0.46265
2.71828
-1.25563
1.85914
-0.396491
2
1.14201
0.320637
2.00115
-0.538504
1.57158
-0.108933
3
1.22571
0.236936
1.79847
-0.335825
1.51209
-0.049445
4
1.27589
0.186756
1.70546
-0.242814
1.49068
-0.028029
5
1.30883
0.153824
1.65248
-0.189833
1.48065
-0.018005
6
1.33199
0.130661
1.61837
-0.155719
1.47518
-0.012529
7
1.34913
0.113518
1.5946
-0.131951
1.47187
-0.009216
8
1.36232
0.10033
1.5771
-0.114455
1.46971
-0.007062
9
1.37277
0.0898766
1.56369
-0.101044
1.46823
-0.005583
10
1.38126
0.0813894
1.55309
-0.090439
1.46717
-0.004525
11
1.38829
0.0743629
1.54449
-0.081844
1.46639
-0.003741
12
1.3942
0.0684508
1.53739
-0.074739
1.46579
-0.003144
13
1.39924
0.0634079
1.53142
-0.068768
1.46533
-0.002680
14
1.40359
0.059056
1.52633
-0.063679
1.46496
-0.002311
15
1.40739
0.0552623
1.52194
-0.059290
1.46466
-0.002014
divisiones
simpson 1/3
error absoluto
1
1.23943
0.22322
2
1.47573
-0.01308
3
1.2564
0.20625
4
1.46371
-0.00106
5
1.31699
0.14566
6
1.46287
-0.00022
7
1.35152
0.11113
8
1.46272
-7E-05
9
1.37309
0.08956
10
1.46268
-3E-05
11
1.38774
0.07491
12
1.46267
-2E-05
13
1.3983
0.06435
14
1.46266
-1E-05
15
1.40626
0.05639

Publicado por en

No hay comentarios:

Publicar un comentario

Suscribirse a: Enviar comentarios (Atom)