0

My problem is as it says in the title, I am trying to use the derivative (with respect to v) of the modified Bessel function of the second kind K_v(x) but with no success.

I read in one of the documentation that besselDK(v,x) would work as a derivative, apparently this is not a recognized function in R. I tried to use the expansion for the derivative, namely 

besselK(v,x)*(1- (1/2v) -log(e*x/2v))

but this doesn't work to give me the correct plot as well. I am trying to plot a function which includes this.

P <- function(x) (1/2)*log(exp(1)/(2*pi*x^(2)))+(3*exp(1/x^(2))/(sqrt(2*pi*x^(2))))*besselK((1/x^(2)),1/2)*(log(exp(1)/x^(2)))
x <- seq(0.1,2,0.01)
plot(x, P(x), xlim=c(0,2), ylim=c(0,1.2), type="l")

From the code above, I get a straight line as a plot. In the correct plot, it should be a curve bending between 1 and 1.5, could someone please tell me the right way to go about it?

Stéphane Laurent
  • 75,186
  • 15
  • 119
  • 225
14thTimeLord
  • 363
  • 1
  • 14
  • `sos::findFn("besselDK")` finds such a function in the [fAsianOptions package](http://finzi.psych.upenn.edu/R/library/fAsianOptions/html/BesselFunctions.html) – Ben Bolker Aug 13 '19 at 19:51
  • Thank you, after loading the package, it did recognize the function 'BesselDK" but it wants the second argument (nu) to be an integer for this function to work, which for me is (1/2). – 14thTimeLord Aug 13 '19 at 20:21
  • This is code taken from the algorithms in Abramowitz and Stegun. You'd have to look to see whether there is something intrinsic about AS's algorithm that is restricted to integer order parameters, or alternatively whether the code could be hacked/tweaked to handle non-integer orders (the documentation does say "an integer value greater or equal to zero" ...) – Ben Bolker Aug 13 '19 at 20:26

1 Answers1

0

The derivative at nu = 1/2 is given here.

f <- function(nu,x){
  besselK(x, nu)
}

library(gsl) # for expint_E1
fprime <- function(x){
  sqrt(pi/2/x) * expint_E1(2*x) * exp(x)
}

nu <- 1/2
h <- 1e-6
x <- 2
(f(nu+h, x) - f(nu,x)) / h 
## [1] 0.02474864
fprime(x)
## [1] 0.02474864
Stéphane Laurent
  • 75,186
  • 15
  • 119
  • 225