Численное интегрирование свертки интерполирующей функции

Question

У меня есть x и y массивы точек данных, которые я использовал для создания интерполяционной функции func_spline как показано ниже.

 import numpy
 from numpy import loadtxt
 from scipy.interpolate import *
 x_given = numpy.arange(1,21400,21400/25000)
 y_given  = loadtxt("Yvalues.txt", comments="#", delimiter=",", unpack=False)

 func_spline = interp1d(x_given,y_given, kind='cubic')
 y_is = func_spline(x_given)

 def alphasinterp(ktsq):
 return func_spline(ktsq)

 import sympy as sp
 import scipy.integrate.quadrature as sciquad

 result = sciquad(alphasinterp,0,21400)
 print(result)

Код успешно выполняет интеграцию, однако я хотел бы изменить код, чтобы разрешить интеграцию вида

   result = sciquad(alphasinterp*f1,0,21400)

где f1 это любая функция (функция ktsq и другие переменные, которые не участвуют в интеграции), что я вьющусь с alphasinterp в интеграции. Например, для конкретного f1 Я получаю ошибку

TypeError: unsupported operand type(s) for *: 'function' and 'FunctionClass'

Как решить? Спасибо! (Как видно из кода, массив y содержит около 21000 точек, поэтому копирование и вставка моих данных здесь, вероятно, недопустимо или нежелательно. Я рад загрузить текстовый файл 'Yvalues.txt', содержащий данные, но я пока не вижу способа сделать это)

Answer 1

Вы должны определить произведение двух функций как другую функцию, которую вы можете использовать при интеграции. Так что в вашем коде это будет выглядеть так:

def alphasinterp(ktsq):
    return func_spline(ktsq)

def f1(ktsq, a1, a2):
    return a1*ktsq+a2# some value

def f_product(ktsq, a1, a2):
    return alphasinterp(ktsq)*f1(ktsq, a1, a2)

def integrated_f(a1, a2):
    return sciquad(f_product,0,21400, args=(a1, a2))

a1=5.0 # just some numbers
a2=3.2
result = integrated_f(a1, a2)

Если вы хотите вычислить свертку, вы должны сделать еще один шаг вперед. С (f*g)(x)=\int f(t)g(xt) dt это было бы что-то вроде

def conv_without_int(t, x):
    return alphasinterp(t)*f1(x-t)
def convolution(x):
    return sciquad(conv_without_int,0,21400, args=(x))

Вы можете сократить это, используя лямбда-нотацию