0

I wrote a function for these equations

f(r) = 1/r(-1/2*((ln(r/rp)/sp)^2)

Y=integral1(f(r)*(1-x/r)^2)/integral2(f(r))

with integral 1 from x to infinity and integral 2 from 0 to infinity

Later I want to fit according to the function but currently I am only using it to create Y for selected x Data.

However, when plotting using certain values for the variables rp and sp (rp: 20.74, sp: 0.053)) the values for Y are completely off from the actual curve for low x values: the values of Y drop to almost 0 for x smaller than approx 0.6.

When changing the variables, e.g. to rp: 20.74 and sp: 0.53 it works.

Graph of the curve with outlier can be found here and the graph of the curve without any outlier (rp: 20.74 and sp: 0.53) can be found here

I assume it has to do with my code respectively the way I am using the integrals, but I can't find out what is the mistake.

import numpy as np
import scipy.integrate
import scipy.optimize
import matplotlib.pyplot as plt

def func(x,rp,sp):
    #Funktion für Kurvenfit
    fn_up = lambda r,rdex : (1/r*np.exp(-1/2*((np.log(r/rp))/sp)**2))*((1-(rdex/r))**2)
    fn_low = lambda r : (1/r*np.exp(-1/2*((np.log(r/rp))/sp)**2))
    num_up = np.asarray([scipy.integrate.quad(fn_up, _x, np.inf,args=(_x,))[0] for _x in x])
    num_low = np.asarray([scipy.integrate.quad(fn_low, 0, np.inf)[0]])
    Kd=num_up/num_low
    return Kd

rp=20.74
sp=0.053
R = np.arange(0.,40.,0.25)
Kd = func(R,rp,sp)
Louise
  • 1
  • 2
    Maybe this isn't a FAQ, but it is at least an OAQ (Occasionally Asked Question). Take a look at my answers to these questions: https://stackoverflow.com/questions/29179778/scipy-quad-uses-only-1-subdivision-and-gives-wrong-result, https://stackoverflow.com/questions/47193045/why-does-integrate-quadlambda-x-xexp-x2-2-sqrt2pi-0-0-100000-give-0, https://stackoverflow.com/questions/37656971/scipy-integrate-pseudo-voigt-function-integral-becomes-0, https://stackoverflow.com/questions/40536467/scipy-integrate-quad-not-returning-the-expected-values – Warren Weckesser Jan 08 '19 at 16:01
  • Thanks. This is plausible and using the points=[] and splitting the integral into the part that contains the signal and the remaining section got my code to work - that is until my sp got smaller again. num_up1 = np.asarray([scipy.integrate.quad(fn_up, _x, 100.,args=(_x,),points=[_x+1.,breakpoint,100.])[0] for _x in x]) with breakpoint = x[-1]+1. and a second num_up2 for integration from 100. to _x I also tried other to integrate to 1000. or 10000. before - with the same result. – Louise Jan 09 '19 at 08:35
  • I suggest plotting the integrand functions for the problematic cases so you can verify that you are breaking up the interval and using `points` appropriately. Additionally (or alternatively), find a formula for the location of the peak of the integrand as a function of the other parameters, and use that in `points` or to break up the integral in a way that ensures that `quad` "sees" the peak. – Warren Weckesser Jan 09 '19 at 17:14
  • thanks, was a bit tricky as my rp value is smaller than my _x but using try/except I got it to work. – Louise Jan 10 '19 at 15:53

0 Answers0