Anonymous user
Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Numerical integration/Gauss-Legendre Quadrature (view source)
Revision as of 18:12, 13 July 2011
, 12 years agoOCaml implementation
(Sorry for breaking the narrative flow) |
(OCaml implementation) |
||
Line 225:
(int #'(lambda (x) x) 5 0 6000)
1.8000000000000004d7</lang>
=={{header|OCaml}}==
<lang OCaml>let rec leg n x = match n with (* Evaluate Legendre polynomial *)
| 0 -> 1.0
| 1 -> x
| k -> let u = 1.0 -. 1.0 /. float k in
(1.0+.u)*.x*.(leg (k-1) x) -. u*.(leg (k-2) x);;
let leg' n x = match n with (* derivative *)
| 0 -> 0.0
| 1 -> 1.0
| _ -> ((leg (n-1) x) -. x*.(leg n x)) *. (float n)/.(1.0-.x*.x);;
let approx_root k n = (* Reversed Francesco Tricomi: 1 <= k <= n *)
let pi = acos (-1.0) and s = float(2*n)
and t = 1.0 +. float(1-4*k)/.float(4*n+2) in
(1.0 -. (float (n-1))/.(s*.s*.s))*.cos(pi*.t);;
let rec refine r n = (* Newton-Raphson *)
let r1 = r -. (leg n r)/.(leg' n r) in
if abs_float (r-.r1) < 2e-16 then r1 else refine r1 n;;
let root k n = refine (approx_root k n) n;;
let node k n = (* Abscissa and weight *)
let r = root k n in
let deriv = leg' n r in
let w = 2.0/.((1.0-.r*.r)*.(deriv*.deriv)) in
(r,w);;
let nodes n =
let rec aux k = if k > n then [] else node k n :: aux (k+1)
in aux 1;;
let quadrature n f a b =
let f1 x = f ((x*.(b-.a) +. a +. b)*.0.5) in
let eval (x,w) = w*.(f1 x) in
0.5*.(b-.a)*.(List.fold_left (+.) 0.0 (List.map eval (nodes n)));;</lang>
An example of how to use it:
<lang OCaml>let calc n =
Printf.printf
"Gauss-Legendre %2d-point quadrature for exp over [-3..3] = %.16f\n"
n (quadrature n exp (-3.0) 3.0);;
calc 5;;
calc 10;;
calc 15;;
calc 20;;</lang>
which outputs
<pre>
Gauss-Legendre 5-point quadrature for exp over [-3..3] = 20.0355777183855608
Gauss-Legendre 10-point quadrature for exp over [-3..3] = 20.0357498548197839
Gauss-Legendre 15-point quadrature for exp over [-3..3] = 20.0357498548198052
Gauss-Legendre 20-point quadrature for exp over [-3..3] = 20.0357498548198052
</pre>
=={{header|J}}==
| |||