1

I have x and y arrays of data points that I have used to create an interpolating function func_spline as shown below.

 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)

The code successfully does the integration however I would like to amend the code to allow integrations of the form

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

where f1 is any function (function of ktsq and other variables that do not partake in integration) that I convolute with alphasinterp in the integration. E.g for a particular f1 I obtain the error

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

How to resolve? Thanks! (as can be seen from the code, the y array contains around 21000 points so a copy and paste of my data here is probably not permissible or advisable. I am happy to upload the text file 'Yvalues.txt' containing the data but I don't see a way to do this yet)

CAF
  • 329
  • 3
  • 14

1 Answers1

1

You have to define the product of the two functions as another function which you can use in the integration. So in your code it would look like:

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)

If you want to compute the convolution you have to take this a step further. With (f*g)(x)=\int f(t)g(x-t) dt it would be something like

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

You could shorten that by using the lambda notation

CodeZero
  • 1,649
  • 11
  • 18
  • Yep I just realised so I deleted my comment but you saw it already :) The upper part of your answer is enough for me so thanks, but my `f1` is in principle a function of the integration variable `ktsq` as well as other variables as noted in my question. I tried to add in other variables but the error is `TypeError: cannot determine truth value of Relational` – CAF Jan 25 '18 at 13:51
  • I think since I am calling a numerical integration method (quad) then I have to assign the additional variables to be numeric rather than symbolically defined but I have not done this before in python (just a thought) – CAF Jan 25 '18 at 13:53
  • Yes, you can only use numeric parameters in the function if you use a numerical intergration method and cannot mix it with symbolic integration or such things. – CodeZero Jan 25 '18 at 13:59
  • Ok great thanks , can you amend your answer to show how this is done in practice? I’ve done it only in mathematica before – CAF Jan 25 '18 at 14:05
  • What do you mean? Extra arguments? It is in the second part of my answer. The x is an extra parameter in the integration – CodeZero Jan 25 '18 at 14:07
  • Yes I just meant like my f1 is now ‘f1(ktsq,a1,a2)’ for example with the last two variables needing to be numerically defined. – CAF Jan 25 '18 at 14:22
  • Ah woops sorry I meant I had to define the additional variables numeric but otherwise left arbitrary. The integration is part of a definition of a chi square function which I’ll minimise next for optimal parameters a1,a2 so I can’t assign them values. If it’s preferrred I’ll post this in another question, thanks. – CAF Jan 25 '18 at 14:33
  • Now you can use the function integrated_f for optimization. – CodeZero Jan 25 '18 at 14:38
  • great thanks, I'm having a little bother using this method in practice for my actual f1, but I put this in another question. – CAF Jan 26 '18 at 12:27