Ламбификация вложенных параметрических интегралов - PullRequest
Новая
       65

Ламбификация вложенных параметрических интегралов

0 голосов
/ 10 сентября 2018

У меня есть функция потерь L (u1, ..., un), которая принимает форму 

L(u) = Integral( (I**2 + J**2) * h, (t,c,d) ) //h=h(t), I=I(t,u)

, где 

I = Integral( f, (x,a,b) )   // f=f(x,t,u)
J = Integral( g, (x,a,b) )   // g=g(x,t,u)

Я хочу численно минимизировать L с помощью scipy, поэтому мне нужно lambdify выражение.

Однако на данный момент lambdify изначально не поддерживает перевод интегралов. С помощью некоторых приемов можно заставить его работать с одиночными параметрическими интегралами, см. Lambdify A Parametric Integral . Однако я не понимаю, как предлагаемое решение могло бы быть распространено на этот обобщенный случай.

Одна идея, которая в принципе должна работать, заключается в следующем:

Возьмем вычислительный граф, определяющий L. Рекурсивно, начиная с листьев, заменяйте каждую символическую операцию соответствующей числовой операцией, выраженной как функция lambda. Однако это приведет к огромному вложению лямбда-функции, которая, как я подозреваю, очень плохо влияет на производительность.

В идеале мы бы хотели получить тот же результат, что и при ручном изготовлении: 

L = lambda u: quad(lambda t: (quad(lambda x: f,a,b)[0]**2 
                  + quad(lambda x: g,a,b)[0]**2)*h, c, d)[0]

MWE: (с использованием кода из старого потока) 

from sympy import *
from scipy.integrate import quad
import numpy as np
def integral_as_quad(function, limits):
    x, a, b = limits
    param = tuple(function.free_symbols - {x})
    f = sp.lambdify((x, *param), function, modules=['numpy'])
    return quad(f, a, b, args=param)[0]

x,a,b,t = sp.symbols('x a b t')
f = (x/t-a)**2
g = (x/t-b)**2
h = exp(-t)
I = Integral( f**2,(x,0,1))
J = Integral( g**2,(x,0,1))
K = Integral( (I**2 + J**2)*h,(t,1,+oo))
F = lambdify( (a,b), K, modules=['sympy',{"Integral": integral_as_quad}])
L = lambda a,b: quad(lambda t: (quad(lambda x: (x/t-a)**2,0,1)[0]**2 
     + quad(lambda x: (x/t-b)**2,0,1)[0]**2)*np.exp(-t), 1, np.inf)[0] 
display(F(1,1))
display(type(F(1,1)))
display(L(1,1))
...