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

У меня есть массивы 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», содержащий данные, но я пока не вижу способа сделать это)


python-3.x interpolation integration cubic-spline numerical-integration
person CAF    schedule 25.01.2018    source источник

Ответы (1)


arrow_upward
1
arrow_downward

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

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(x-t) 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))

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

person datasailor    schedule 25.01.2018
comment
Да, я только что понял, поэтому удалил свой комментарий, но вы его уже видели :) Верхней части вашего ответа мне достаточно, так что спасибо, но мой f1 в принципе является функцией переменной интеграции ktsq, а также других переменных, как указано в мой вопрос. Я пытался добавить другие переменные, но ошибка TypeError: cannot determine truth value of Relational - person CAF; 25.01.2018
comment
Я думаю, что, поскольку я вызываю метод численного интегрирования (quad), мне нужно назначить дополнительные переменные числовыми, а не символическими, но я не делал этого раньше в python (просто мысль) - person CAF; 25.01.2018
comment
Да, вы можете использовать числовые параметры в функции только в том случае, если вы используете метод числового интегрирования и не можете смешивать его с символьным интегрированием или подобными вещами. - person datasailor; 25.01.2018
comment
Хорошо, большое спасибо, можете ли вы изменить свой ответ, чтобы показать, как это делается на практике? Раньше я делал это только в математике - person CAF; 25.01.2018
comment
Что ты имеешь в виду? Лишние аргументы? Это во второй части моего ответа. x является дополнительным параметром в интегрировании - person datasailor; 25.01.2018
comment
Да, я просто имел в виду, что моя f1 теперь «f1 (ktsq, a1, a2)», например, с двумя последними переменными, которые необходимо определить численно. - person CAF; 25.01.2018
comment
Ах, извините, я имел в виду, что должен был определить дополнительные числовые переменные, но в остальном оставил произвольные. Интегрирование является частью определения функции хи-квадрат, которую я минимизирую далее для оптимальных параметров a1,a2, поэтому я не могу присвоить им значения. Если это предпочтительнее, я опубликую это в другом вопросе, спасибо. - person CAF; 25.01.2018
comment
Теперь вы можете использовать функцию Integrated_f для оптимизации. - person datasailor; 25.01.2018
comment
Большое спасибо, у меня есть небольшие затруднения с использованием этого метода на практике для моей реальной f1, но я поставил это в другой вопрос. - person CAF; 26.01.2018