Programación de macros VBA para Excel M.Sc. Walter Mora F., M.Sc. José Luis Espinoza B.

Integral doble gaussiana.

Figure 47: Integral doble gaussiana.

Si queremos aplicar la cuadratura gaussiana a $\,\displaystyle{\int_{a}^{b}\int_{c(x)}^{d(x)}f(x,y)\,dy\,dx}\,$ se debe hacer, primero que todo, un cambio de variable para cambiar $\,[c(x), \, d(x)]\,$ a $\,[-1,1]\,$ . Esto se hace de manera parecida al procedimiento numérico.

$\,\displaystyle{\int_{a}^{b}\int_{c(x)}^{d(x)}f(x,y)\,dy\,dx} \approx \displays... ...}c_{nj}f\left[ x, \, 0.5\cdot(d(x)-c(x)\,x_{nj}\,+\,d(x)+c(x)) \, dx \right]}\,$

y luego debemos volver hacer esto para $\,[a,b]\,$ .

Cada Integral puede ser aproximada usando las raíces de un polinomio de Legendre de grado distinto.

El algoritmo es el siguiente

Entrada: $\,a,b\,$ y los enteros $\,m,n\,$ correspondientes al grado del polinomio de Legendre que se usa en cada integral. Deben ser enteros adecuados de acuerdo a los valores $\,c_{ij},\,x_{ij}\,$ disponibles.
Tome

$\,h_1=(b-a)/2\,$ ;

$\,h_2 =(b+a)/2\,$

$\,J=0\,$
Para
hacer
1. Tome $\,JX = 0\,$
  
  $\,x=h_1x_{mi}+h_2\,$
  
  $\,d_1 =d(x)\,$
  
  $\,c_1=c(x)\,$
  
  $\,k_1=(d_1-c_1)/2\,$
  
  $\,k_2= (d_1+c_1)/2\,$
2. para $\,j=1,2,...,n \;\,$ hacer
  
  $\,y = k_1x_{nj}+k_2\,$
  
  $\,Q=f(x,y)\,$
  
  $\,JX=JX+c_{nj}Q\,$
3. Tome $\,J=J+c_{mi}k_1JX\,$
Tome $\,J = h_1 J\,$
Salida: la aproximación de la integral es $\,J\,$

En la implementación debemos evaluar funciones de una variable, $\,c(x)\,$ y $\,d(x)\,$ ,y una función de dos variables, $\,Q=f(x,y)\,$ . Para hacer esto, almacenamos la expresión actual y procedemos a evaluar ya sea con Eval1() para funciones de una variable o Eval() para varias variables. En el código aparecen las siguientes líneas para este fin

OK = Fun.StoreExpression(dx)      'formula actual es 'd(x)'
d1 = Fun.Eval1(xa) 'd(xa)
    ...
OK = Fun.StoreExpression(cx)      'formula actual es 'c(x)'
c1 = Fun.Eval1(xa) 'c(xa)
    ...
OK = Fun.StoreExpression(fxy)      'formula actual es 'f(x,y)'
Fun.Variable("x") = xa
Fun.Variable("y") = ya
Q = Fun.Eval() 'retorna f(xa,ya)
    ...

El código completo del botón es (faltarían xi( ) y yi( ) )

Private Sub CommandButton1_Click()


Dim a, b, h1, h2, Jt, JX, xa, d1, c1, k1, k2, ya, Q As Double
Dim n, m As Integer
Dim fxy, cx, dx As String 'funciones
Dim OK As Boolean
Dim Fun As New clsMathParser      ' así se llama el módulo de clase aquí

fxy = Cells(2, 2)
cx = Cells(3, 5)
dx = Cells(3, 6)
a = Cells(3, 3)
b = Cells(3, 4)
n = Cells(3, 7)
m = Cells(3, 8)

Cells(7, 9) = "" 'limpiar cálculos anteriores

If n < 2 Or n > 5 Or m < 2 Or m > 5 Then
    MsgBox ("Solo tenemos datos para n y m entre 2 y 5")
    Exit Sub
End If
'paso 1
h1 = (b - a) / 2
h2 = (b + a) / 2
Jt = 0
'paso2
For i = 1 To m
    'paso 3
    JX = 0
    xa = h1 * xi(m, i) + h2
    OK = Fun.StoreExpression(dx)      'formula actual es 'd(x)'
    If Not OK Then GoTo Error_Handler

    d1 = Fun.Eval1(xa) 'd(xa)

    OK = Fun.StoreExpression(cx)      'formula actual es 'c(x)'
    If Not OK Then GoTo Error_Handler

    c1 = Fun.Eval1(xa) 'c(xa)

    k1 = (d1 - c1) / 2
    k2 = (d1 + c1) / 2

    For J = 1 To n
        ya = k1 * xi(n, J) + k2

        OK = Fun.StoreExpression(fxy)      'formula actual es 'f(x,y)'
        If Not OK Then GoTo Error_Handler
        Fun.Variable("x") = xa
        Fun.Variable("y") = ya
        Q = Fun.Eval() 'retorna f(xa,ya)

        JX = JX + ci(n, J) * Q
    Next J

    Jt = Jt + ci(m, i) * k1 * JX

Next i

Jt = h1 * Jt

Cells(7, 9) = Jt  ' aproxx de la integral

'---------------------------------------------------------------
If Err Then GoTo Error_Handler
Error_Handler: Debug.Print Fun.ErrorDescription
'---------------------------------------------------------------

End Sub

Nota: Observe que como $\,\displaystyle{\int _{a}^{b}\int _{0}^{1}\,2\,y\,f(x)\,dy\,dx=\int _{a}^{b}f(x)\,dx},\,$ entonces esta implementación es adecuada para la cuadratura gaussiana en una variable.

Nota: Se puede aplicar la misma idea y aproximar una integral doble usando Simpson.