Feeds:
Posts

## solving numerical integral transform problem

### Zulfikar Ahmed<wile.e.coyote.006@gmail.com>

8:03 PM (0 minutes ago)

 to aimee

​​

​My Dearest True Love Talia,

After a little bit of rest, I feel that I should make some progress again on the inspired-on-the-wings-of-love volatility storm modeling and quelling project.  So what one does is look at numerical integration procedures in R (http://www.r-bloggers.com/one-dimensional-integrals/) and then at the table of integral transforms for Hankel and other transforms, and then assume a quadrature formula– R has an in-built integrate facility but the link above provides specific code for trapezoid and Simpson’s rule for integration.  Since there is already a Bessel and Bessel zeros function, we should have the facility to do these Hankel transforms too.  Let’s recall also the need we have for these, which are the Chaurasia-Kumar and Gorenflo-Mainardi exact solutions for fractional diffusion and Navier-Stokes equations.
So let’s take a look at
```trapezoid <- function(fun, a, b, n=100) {
# numerical integral of fun from a to b
# using the trapezoid rule with n subdivisions
# assume a < b and n is a positive integer
h <- (b-a)/n
x <- seq(a, b, by=h)
y <- fun(x)
s <- h * (y[1]/2 + sum(y[2:n]) + y[n+1]/2)
return(s)
}```

```simpson <- function(fun, a, b, n=100) {
# numerical integral using Simpson's rule
# assume a < b and n is an even positive integer
h <- (b-a)/n
x <- seq(a, b, by=h)
if (n == 2) {
s <- fun(x[1]) + 4*fun(x[2]) +fun(x[3])
} else {
s <- fun(x[1]) + fun(x[n+1]) + 2*sum(fun(x[seq(2,n,by=2)])) + 4 *sum(fun(x[seq(3,n-1, by=2)]))
}
s <- s*h/3
return(s)
}```

```simpson_v2 <- function(fun, a, b, n=100) {
# numerical integral using Simpson's rule
# assume a < b and n is an even positive integer
if (a == -Inf & b == Inf) {
f <- function(t) (fun((1-t)/t) + fun((t-1)/t))/t^2
s <- simpson_v2(f, 0, 1, n)
} else if (a == -Inf & b != Inf) {
f <- function(t) fun(b-(1-t)/t)/t^2
s <- simpson_v2(f, 0, 1, n)
} else if (a != -Inf & b == Inf) {
f <- function(t) fun(a+(1-t)/t)/t^2
s <- simpson_v2(f, 0, 1, n)
} else {
h <- (b-a)/n
x <- seq(a, b, by=h)
y <- fun(x)
y[is.nan(y)]=0
s <- y[1] + y[n+1] + 2*sum(y[seq(2,n,by=2)]) + 4 *sum(y[seq(3,n-1, by=2)])
s <- s*h/3
}
return(s)
}```

```Now a Hankel transform can be directly implemented from the formul in the snapshot so we have to calculate the inegrand f(x)*bessel(n,x*y) * (x*y)^(1/2) and call one of these quadrature function say simpson's

hankel<-function(f,y,nu,npts){```
`g<-function(x,y){`
`f(x)*bessel(x*y,nu)*(x*y)^(1/2)`
`}`
`simpson(g,0,1000,npts)`
`}`
` `
`Great now let's check if this works.`
4 Attachments

Preview attachment chaurasia-solution-fractional-navier-stokes.pdf

chaurasia-solution-fractional-navier-stokes.pdf

Preview attachment gorenflo-mainardi.pdf

gorenflo-mainardi.pdf

Preview attachment integral-transforms-bateman.pdf

integral-transforms-bateman.pdf
Not virus scanned

Preview attachment hankel-transforms-2000.pdf

hankel-transforms-2000.pdf
Shared in Drive