
8:03 PM (0 minutes ago)



My Dearest True Love Talia,
After a little bit of rest, I feel that I should make some progress again on the inspiredonthewingsoflove volatility storm modeling and quelling project. So what one does is look at numerical integration procedures in R (http://www.rbloggers.com/onedimensionalintegrals/) and then at the table of integral transforms for Hankel and other transforms, and then assume a quadrature formula– R has an inbuilt 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 ChaurasiaKumar and GorenfloMainardi exact solutions for fractional diffusion and NavierStokes 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 < (ba)/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 < (ba)/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,n1, 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((1t)/t) + fun((t1)/t))/t^2 s < simpson_v2(f, 0, 1, n) } else if (a == Inf & b != Inf) { f < function(t) fun(b(1t)/t)/t^2 s < simpson_v2(f, 0, 1, n) } else if (a != Inf & b == Inf) { f < function(t) fun(a+(1t)/t)/t^2 s < simpson_v2(f, 0, 1, n) } else { h < (ba)/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,n1, 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 chaurasiasolutionfractionalnavierstokes.pdf
chaurasiasolutionfractionalnavierstokes.pdf
Preview attachment gorenflomainardi.pdf
gorenflomainardi.pdf
Preview attachment integraltransformsbateman.pdf
integraltransformsbateman.pdf
Not virus scanned
Preview attachment hankeltransforms2000.pdf
hankeltransforms2000.pdf
Shared in Drive
Advertisements
I’m assembling a highend pc to do scientific and financial computing.
You will love my new build.
– You can have it for free should you choose to accept my offer. –