Tengo el siguiente problema. Tengo una función f definida en python usando funciones numpy. La función es suave e integrable en reales positivos. Quiero construir la antiderivada doble de la función (asumiendo que tanto el valor como la pendiente de la antiderivada en 0 son 0) para poder evaluarla en cualquier real positivo menor que 100.
Definición de antiderivada de f en x :
integrate f(s) with s from 0 to x Definición de antiderivada doble de f en x :
integrate (integrate f(t) with t from 0 to s) with s from 0 to x La forma real de f no es importante, así que usaré una simple por conveniencia. Pero tenga en cuenta que aunque mi ejemplo tiene una forma cerrada conocida, mi función real no la tiene.
import numpy as np f = lambda x: np.exp(-x)*xMi solución es construir la antiderivada como una matriz usando integración numérica ingenua:
N = 10000 delta = 100/N xs = np.linspace(0,100,N+1) vs = f(xs) avs = np.cumsum(vs)*delta aavs = np.cumsum(avs)*delta Esto, por supuesto, funciona, pero me da matrices en lugar de funciones. Pero esto no es un gran problema ya que puedo interpolar aavs usando una spline para obtener una función y deshacerme de las matrices.
from scipy.interpolate import UnivariateSpline aaf = UnivariateSpline(xs, aavs) La función aaf es aproximadamente la doble antiderivada de f .
El problema es que aunque funciona, hay un poco de sobrecarga antes de que pueda obtener mi función y la precisión es costosa.
Mi otra idea era interpolar f por un spline y tomar la antiderivada de eso, sin embargo, esto introduce errores numéricos que son demasiado grandes para lo que quiero usar en la función.
¿Hay alguna manera mejor de hacer eso? Por mejor quiero decir más rápido sin sacrificar la precisión.
Editar: lo que espero que sea posible es usar algún tipo de transformada de Fourier para evitar la integración dos veces. Espero que haya alguna transformación conveniente de vs que permita multiplicar los valores por componentes con xs y transformar de nuevo para obtener la antiderivada doble. Jugué un poco con esto, pero me perdí.
Editar: descubrí que al usar la regla trapezoidal en lugar de una suma ingenua, aumenta bastante la precisión. El uso de la regla de Simpson debería aumentar aún más la precisión, pero es un poco complicado hacerlo con matrices numpy.
Editar: como @ user202729 se queja legítimamente, esto parece estar mal. La razón por la que parece mal es porque me he saltado algunos detalles. Explico aquí por qué lo que digo tiene sentido, pero no afecta mi pregunta.
Mi objetivo real no es encontrar la antiderivada doble de f , sino encontrar una transformación de esto. Me he saltado eso porque creo que solo confunde el asunto.
La función f decae exponencialmente cuando x se acerca a 0 o infinito. Estoy minimizando el error numérico en la integración comenzando la suma desde 0 y subiendo hasta aproximadamente el pico de f . Esto asegura que el error relativo sea aproximadamente constante. Luego empiezo desde la dirección opuesta desde una x muy grande y vuelvo al pico. Luego hago lo mismo con los valores de las antiderivadas.
Luego transformo el aavs por otra función que es sensible a los errores numéricos. Luego encuentro la región donde los errores son grandes (los valores oscilan violentamente) y dejo caer estos valores. Finalmente, aproximo lo que creo que son buenos valores mediante una spline.
Ahora, si uso spline para aproximar f , introduce un error absoluto que es el término dominante en un intervalo bastante grande. Esto se "integra" dos veces y termina siendo un error relativo bastante grande en aavs . Luego, una vez que transformo aavs , encuentro que la 'buena región' se ha reducido considerablemente.
EDITAR: La forma real de f es algo que todavía estoy investigando. Sin embargo, va a ser una generalización de la distribución lognormal. Ahora mismo estoy jugando con la siguiente familia.
Comienzo definiendo una generalización de la distribución normal:
def pdf_n(params, center=0.0, slope=8): scale, min, diff = params if diff > 0: r = min l = min + diff else: r = min - diff l = min def retfun(m): x = (m - center)/scale E = special.expit(slope*x)*(r - l) + l return np.exp( -np.power(1 + x*x, E)/2 ) return np.vectorize(retfun)Puede que no sea obvio lo que está sucediendo aquí, pero el resultado es bastante simple. La función decae como exp(-x^(2l)) a la izquierda y como exp(-x^(2r)) a la derecha. Para min=1 y diff=0, esta es la distribución normal. Tenga en cuenta que esto no está normalizado. Entonces defino
g = pdf(params) f = np.vectorize(lambda x:g(np.log(x))/x/area) donde area es la constante de normalización.
Tenga en cuenta que este no es el código real que uso. Lo desnudé al mínimo.
Primero, creé una versión del código que encontré más intuitiva. Aquí multiplico los valores de la suma acumulada por los anchos de los contenedores. Creo que hay un pequeño error en la versión original del código relacionado con el problema del ancho del contenedor.
import numpy as np f = lambda x: np.exp(-x)*x N = 1000 xs = np.linspace(0,100,N+1) domainwidth = ( np.max(xs) - np.min(xs) ) binwidth = domainwidth / N vs = f(xs) avs = np.cumsum(vs)*binwidth aavs = np.cumsum(avs)*binwidthA continuación, para la visualización, aquí hay un código de trazado muy simple:
import matplotlib import matplotlib.pyplot as plt plt.figure() plt.scatter( xs, vs ) plt.figure() plt.scatter( xs, avs ) plt.figure() plt.scatter( xs, aavs ) plt.show()La primera integral coincide con el resultado conocido de la expresión de ejemplo y se puede ver en wolframio
A continuación se muestra una función simple que extrae un elemento de la segunda derivada. Tenga en cuenta que int es una mala función de redondeo. Supongo que esto es lo que ya has implementado.
def extract_double_antideriv_value(x): return aavs[int(x/binwidth)] singleresult = extract_double_antideriv_value(50.24) print('singleresult', singleresult)Cualesquiera que sean los pasos de cálculo completos que se requieran, debemos conocerlos antes de que podamos comenzar a optimizar. ¿Tienes un millón de funciones diferentes para integrar? Si solo necesita consultar una única antiderivada doble muchas veces, su solución original debería ser bastante ideal.
Aproximación simbólica: ¿Ha considerado aproximaciones a la función original f , que puede tener soluciones de integración de forma cerrada? Tiene un dominio limitado en el que vive la función. ¿Quizás aproximar f con una serie de Taylor (que se puede construir con un error máximo conocido ) y luego integrar exactamente? (considere Pade, Taylor, Fourier, Cheby, Lagrange (como lo sugiere otra respuesta), etc.)
Trucos de registro: otra alternativa para lidiar con errores puntiagudos sería tomar el registro de su función original. ¿ f siempre es positivo? ¿El error de integración se debe a que la vecindad alrededor del máximo es muy pequeña? Si es así, puedes estudiar ln(f) o incluso ln(ln(f)) en su lugar. Realmente ayudaría a entender cómo se ve f más.
Trucos de integración de aproximación Existen innumerables trucos de integración en general, que pueden hacer soluciones aproximadas de forma cerrada para integrales que se pueden deshacer. Una muy común cuando se trata de funciones exponenciales (¿creo que la tuya es exponencial?) es usar el método de Laplace . Pero qué truco sacar de la bolsa depende en gran medida de las condiciones que satisfaga f .