trapecio - Escribir una suma infinita(de una función que tiene una integral) en Python

integral trapecio python

Question

Estoy obligado a mostrar que:

Lo molesto es que c_i es igual a la integral de la función G. Aquí está mi intento.

import numpy as np
from scipy.integrate import quad

def G(x,n):
    P = (np.sqrt(735))*(np.sqrt(2))*np.sin(n*np.pi*x)*((x**3.0) - (11.0/7.0)*(x**2.0) + (4.0/7.0)*(x))
    return P

def Sum(x, n):
    i = 1
    S = 0
    I, err = quad(G, 0, 1, args=(i))
    while (i<n):
        S = S + I
        i = i + 1
    return S
x = np.linspace(0, 1, 250)    
print Sum(x, 200)

El problema que estoy teniendo es escribir la parte de la suma. Cuando ejecuto esto obtengo un número mayor cuantos más valores le doy. Si se elige que n es bastante alto (en lugar de infinito), se puede mostrar cómo la suma tiende a 1

Answer 1

Tienes errror en tu resumen. Mueva quad dentro de loop, y evalúe la suma de cuadrados (en lugar de la suma evaluada actualmente (c_i)) de esta manera:

def Sum(x, n):
    S = 0
    for i in range(n):
        I, _err = quad(G, 0, 1, args=(i))
        S = S + I**2
    return S

Por cierto, ¿cómo puedes DEMOSTRAR que el límite de esta suma es igual a 1? Solo se puede mostrar mediante algunas evaluaciones matemáticas, no computaciones. Puedes simplemente ilustrar eso.

Answer 2

Este problema tiene mucho valor pedagógico.

Como señala @alko, este problema se puede resolver analíticamente. Si la intención es demostrar que la suma es igual a uno, entonces debe hacerse analíticamente.

Tal vez, esta es simplemente una versión más simple de algo que debe hacerse y el problema real no se puede resolver analíticamente. En ese caso, resolver un problema más simple como este es un buen primer paso. Desafortunadamente, cada vez que resolvemos algo numéricamente, presentamos nuevos conjuntos de problemas.

Procedamos como se sugiere en la pregunta y la corrección dada por @alko.

import numpy as np
import scipy.integrate as integ

def g(x) :
    return np.sqrt(1470) * (x**3-11./7*x**2+4./7*x)

def G(x,n) :
    return np.sin(n*np.pi*x) * g(x)

def do_integral (N) :
    res = np.zeros(N)
    err = np.zeros(N)
    for n in xrange(1,N) :
        (res[n],err[n]) = integ.quad(G, 0, 1, args=(n,))
    return (res,err)

(c,err) = do_integral(500)
S = np.cumsum(c[1:]**2) # skip n=0

(La razón por la que defino dos funciones "G" será evidente a continuación).

Una vez ejecutado, la matriz S contendrá la suma deseada como una función de n. Debería ser uno. Ahora, cuando ejecutamos esto dentro de ipython, será algo lento y la primera vez recibiremos un mensaje de advertencia largo, pero parece que el código se ejecutará. Además, si lo ejecutamos nuevamente (en la misma sesión de ipython) no recibiremos un mensaje de advertencia por lo que podemos ignorarlo, ¿verdad? Incorrecto pero lo ignoraremos de todos modos, ya que es común hacerlo.

Si observamos S , parece que está mostrando lo que queremos hasta aproximadamente S[200] donde las cosas comienzan a salir mal, ¡el valor comienza a crecer! ¿Qué pasa con la computadora? Nada, e ignoramos otro indicador de un problema. quad() devuelve una estimación de error junto con una estimación de la integral. Por lo general, ignoramos el error estimado, pero no deberíamos. Si graficamos los valores de S y error, encontramos lo siguiente.

Entonces vemos que sí, el valor de S sale terriblemente mal, pero quad() también nos decía que esto sucedería. De hecho, la advertencia que ignorábamos también nos decía lo mismo.

¿Cómo lo entendemos y lo solucionamos? En este punto, me detendré con la narración de historias. Si mirar fijamente el G(x,n) no deja claro, sería bueno trazar esa función para un n grande sobre el rango de integración. Descubriremos que es una función tremendamente oscilante, por lo que no es sorprendente que sea difícil de integrar numéricamente. Debe haber una mejor manera.

Por supuesto que hay una mejor manera. Si miramos la documentación para quad() y ejecutamos quad_explain() podemos aprender sobre las funciones de peso. El seno es una función de peso común que aparece en integrales por lo que existen técnicas especiales para manejar este caso. Por lo tanto, un mejor enfoque es el siguiente (ahora vemos por qué definí g(x) :

def do_integral_weighted (N) :
    res = np.zeros(N)
    err = np.zeros(N)
    for n in xrange(1,N) :
        (res[n],err[n]) = integ.quad(g, 0, 1, weight=''sin'', wvar=n*np.pi)
    return (res,err)

(cw,errw) = do_integral_weighted(500)
Sw = np.cumsum(cw[1:]**2) # skip n=0

Esto nos parece que se ejecuta mucho más rápido y es mucho más preciso como se muestra en la trama Entonces obtenemos una respuesta estable, sin advertencias impresas y un pequeño error de integración.

Aprendemos algunas cosas al trabajar este problema

1. Los problemas que pueden resolverse analíticamente deben resolverse analíticamente.
2. La implementación numérica de una expresión matemática a menudo no es tan directa como podríamos haber esperado.
3. No ignore las advertencias y no ignore la información de error proporcionada por las rutinas que estamos utilizando.
4. Las técnicas numéricas son diferentes a las técnicas analíticas. numpy / scipy proporciona herramientas extremadamente potentes. Necesitamos explorar estas herramientas para explotar completamente su poder.