Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Numerical integration/Gauss-Legendre Quadrature (view source)
Revision as of 15:07, 10 April 2024
, 1 month agoAdded FreeBASIC
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
(Added FreeBASIC) |
||
| (5 intermediate revisions by 3 users not shown) | |||
Line 128:
integral: 20.035577718
</pre>
=={{header|ATS}}==
{{trans|Common Lisp}}
This is a very close translation of the Common Lisp.
(A lot of the "ATS-ism" is completely optional. For instance, you can use <code>arrszref</code> instead of <code>arrayref</code>, if you want bounds checking at runtime instead of compile-time. But then debugging and regression-prevention become harder, and in that particular case the code will almost surely be slower.
And, if I may grumble a bit: ''Some'' of us ''do not'' think "turning off bounds checking for production" is acceptable. It is at best something to tolerate grudgingly.)
<syntaxhighlight lang="ats">
#include "share/atspre_staload.hats"
%{^
#include <float.h>
#include <math.h>
%}
extern fn {tk : tkind} g0float_pi : () -<> g0float tk
extern fn {tk : tkind} g0float_cos : g0float tk -<> g0float tk
extern fn {tk : tkind} g0float_exp : g0float tk -<> g0float tk
implement g0float_pi<dblknd> () = $extval (double, "M_PI")
implement g0float_cos<dblknd> x = $extfcall (double, "cos", x)
implement g0float_exp<dblknd> x = $extfcall (double, "exp", x)
macdef PI = g0float_pi ()
overload cos with g0float_cos
overload exp with g0float_exp
macdef NAN = g0f2f ($extval (float, "NAN"))
macdef Zero = g0i2f 0
macdef One = g0i2f 1
macdef Two = g0i2f 2
(* Computes the initial guess for the root i of a n-order Legendre
polynomial. *)
fn {tk : tkind}
guess {n, i : int | 1 <= i; i <= n}
(n : int n, i : int i) :<> g0float tk =
cos (PI * ((g0i2f i - g0f2f 0.25) / (g0i2f n + g0f2f 0.5)))
(* Computes and evaluates the degree-n Legendre polynomial at the
point x. *)
fn {tk : tkind}
legpoly {n : pos}
(n : int n, x : g0float tk) :<> g0float tk =
let
fun
loop {i : int | 2 <= i; i <= n + 1} .<n + 1 - i>.
(i : int i, pa : g0float tk, pb : g0float tk)
:<> g0float tk =
if i = succ n then
pb
else
let
val iflt = (g0i2f i) : g0float tk
val pn = (((iflt + iflt - One) / iflt) * x * pb)
- (((iflt - One) / iflt) * pa)
in
loop (succ i, pb, pn)
end
in
if n = 0 then
One
else if n = 1 then
x
else
loop (2, One, x)
end
(* Computes and evaluates the derivative of an n-order Legendre
polynomial at point x. *)
fn {tk : tkind}
legdiff {n : int | 2 <= n}
(n : int n, x : g0float tk) :<> g0float tk =
(g0i2f n / ((x * x) - One))
* ((x * legpoly<tk> (n, x)) - legpoly<tk> (pred n, x))
(* Computes the n nodes for an n-point quadrature rule (the n roots of
a degree-n polynomial). *)
fn {tk : tkind}
nodes {n : int | 2 <= n}
(n : int n) :<!refwrt> arrayref (g0float tk, n) =
let
val x = arrayref_make_elt<g0float tk> (i2sz n, Zero)
fn
v_update (v : g0float tk) :<> g0float tk =
v - (legpoly<tk> (n, v) / legdiff<tk> (n, v))
var i : Int
in
for* {i : nat | i <= n} .<n - i>.
(i : int i) =>
(i := 0; i <> n; i := succ i)
let
val v = guess<tk> (n, succ i)
val v = v_update v
val v = v_update v
val v = v_update v
val v = v_update v
val v = v_update v
in
x[i] := v
end;
x
end
(* Computes the weight for an degree-n polynomial at the node x. *)
fn {tk : tkind}
legwts {n : int | 2 <= n}
(n : int n, x : g0float tk) :<> g0float tk =
(* Here I am having slightly excessive fun with notation: *)
Two / ((One - (x * x)) * (y * y where {val y = legdiff<tk> (n, x)}))
(* Normally I would not write code in such fashion. :) Nevertheless,
it is interesting that this works. *)
(* Takes an array of nodes x and computes an array of corresponding
weights w. Note that x is an arrayref, not an arrszref, and so
(unlike in the Common Lisp) we have to tell the function the size
of the new array w. That information is not otherwise stored AT
RUNTIME. The ATS compiler, however, will force us AT COMPILE TIME
to pass the correct size. *)
fn {tk : tkind}
weights {n : int | 2 <= n}
(n : int n, x : arrayref (g0float tk, n))
:<!refwrt> arrayref (g0float tk, n) =
let
val w = arrayref_make_elt<g0float tk> (i2sz n, Zero)
var i : Int
in
for* {i : nat | i <= n} .<n - i>.
(i : int i) =>
(i := 0; i <> n; i := succ i)
w[i] := legwts (n, x[i]);
w
end
(* Estimates the definite integral of a function on [a,b], using an
n-point Gauss-Legendre quadrature rule. *)
fn {tk : tkind}
quad {n : int | 2 <= n}
(f : g0float tk -<> g0float tk,
n : int n,
a : g0float tk,
b : g0float tk) :<> g0float tk =
let
val x = $effmask_ref ($effmask_wrt (nodes<tk> n))
val w = $effmask_ref ($effmask_wrt (weights<tk> (n, x)))
val ahalf = g0f2f 0.5 * a and bhalf = g0f2f 0.5 * b
val C1 = bhalf - ahalf and C2 = ahalf + bhalf
fun
loop {i : nat | i <= n} .<n - i>.
(i : int i, sum : g0float tk) :<> g0float tk =
if i = n then
sum
else
let
val y = $effmask_ref (w[i] * f ((C1 * x[i]) + C2))
in
loop (succ i, sum + y)
end
in
C1 * loop (0, Zero)
end
implement
main () =
let
val outf = stdout_ref
in
fprintln! (outf, "nodes<dblknd> 5");
fprint_arrayref_sep<double> (outf, nodes<dblknd> (5),
i2sz 5, " ");
fprintln! (outf); fprintln! (outf);
fprintln! (outf, "weights (nodes<dblknd> 5)");
fprint_arrayref_sep<double> (outf, weights (5, nodes<dblknd> (5)),
i2sz 5, " ");
fprintln! (outf); fprintln! (outf);
fprintln! (outf, "quad (lam x => exp x, 5, ~3.0, 3.0) = ",
quad (lam x => exp x, 5, ~3.0, 3.0));
fprintln! (outf);
fprintln! (outf, "More examples, borrowed from the Common Lisp:");
fprintln! (outf, "quad (lam x => x ** 3, 5, 0.0, 1.0) = ",
quad (lam x => x ** 3, 5, 0.0, 1.0));
fprintln! (outf, "quad (lam x => 1.0 / x, 5, 1.0, 100.0) = ",
quad (lam x => 1.0 / x, 5, 1.0, 100.0));
fprintln! (outf, "quad (lam x => x, 5, 0.0, 5000.0) = ",
quad (lam x => x, 5, 0.0, 5000.0));
fprintln! (outf, "quad (lam x => x, 5, 0.0, 6000.0) = ",
quad (lam x => x, 5, 0.0, 6000.0));
0
end
</syntaxhighlight>
{{out}}
<pre>$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW gauss_legendre_task.dats -lgc -lm && ./a.out
nodes<dblknd> 5
0.906180 0.538469 0.000000 -0.538469 -0.906180
weights (nodes<dblknd> 5)
0.236927 0.478629 0.568889 0.478629 0.236927
quad (lam x => exp x, 5, ~3.0, 3.0) = 20.035578
More examples, borrowed from the Common Lisp:
quad (lam x => x ** 3, 5, 0.0, 1.0) = 0.250000
quad (lam x => 1.0 / x, 5, 1.0, 100.0) = 4.059148
quad (lam x => x, 5, 0.0, 5000.0) = 12500000.000000
quad (lam x => x, 5, 0.0, 6000.0) = 18000000.000000</pre>
=={{header|Axiom}}==
Line 857 ⟶ 1,067:
20 20.0357498548198037979491872388495 -.82E-28
</pre>
=={{header|FreeBASIC}}==
{{trans|Wren}}
<syntaxhighlight lang="vbnet">#define PI 4 * Atn(1)
Const As Double LIM = 5
Dim Shared As Double lroots(LIM - 1)
Dim Shared As Double weight(LIM - 1)
Dim Shared As Double lcoef(LIM, LIM)
For i As Integer = 0 To LIM
For j As Integer = 0 To LIM
lcoef(i, j) = 0
Next j
Next i
Sub legeCoef()
lcoef(0, 0) = 1
lcoef(1, 1) = 1
For n As Integer = 2 To LIM
lcoef(n, 0) = -(n - 1) * lcoef(n - 2, 0) / n
For i As Integer = 1 To n
lcoef(n, i) = ((2 * n - 1) * lcoef(n - 1, i - 1) - (n - 1) * lcoef(n - 2, i)) / n
Next i
Next n
End Sub
Function legeEval(n As Integer, x As Double) As Double
Dim As Double s = lcoef(n, n)
For i As Integer = n To 1 Step -1
s = s * x + lcoef(n, i - 1)
Next i
Return s
End Function
Function legeDiff(n As Integer, x As Double) As Double
Return n * (x * legeEval(n, x) - legeEval(n - 1, x)) / (x * x - 1)
End Function
Sub legeRoots()
Dim As Double x = 0
Dim As Double x1 = 0
For i As Integer = 1 To LIM
x = Cos(PI * (i - 0.25) / (LIM + 0.5))
Do
x1 = x
x = x - legeEval(LIM, x) / legeDiff(LIM, x)
Loop Until x = x1
lroots(i - 1) = x
x1 = legeDiff(LIM, x)
weight(i - 1) = 2 / ((1 - x * x) * x1 * x1)
Next i
End Sub
Function legeIntegrate(f As Function (As Double) As Double, a As Double, b As Double) As Double
Dim As Double c1 = (b - a) / 2
Dim As Double c2 = (b + a) / 2
Dim As Double sum = 0
For i As Integer = 0 To LIM - 1
sum = sum + weight(i) * f(c1 * lroots(i) + c2)
Next i
Return c1 * sum
End Function
legeCoef()
legeRoots()
Print "Roots: ";
For i As Integer = 0 To LIM - 1
Print Using " ##.######"; lroots(i);
Next i
Print
Print "Weight:";
For i As Integer = 0 To LIM - 1
Print Using " ##.######"; weight(i);
Next i
Print
Function f(x As Double) As Double
Return Exp(x)
End Function
Dim As Double actual = Exp(3) - Exp(-3)
Print Using !"Integrating exp(x) over [-3, 3]:\n\t########.######,\ncompared to actual\n\t########.######"; legeIntegrate(@f, -3, 3); actual
Sleep</syntaxhighlight>
{{out}}
<pre>Roots: 0.906180 0.538469 0.000000 -0.538469 -0.906180
Weight: 0.236927 0.478629 0.568889 0.478629 0.236927
Integrating exp(x) over [-3, 3]:
20.035578,
compared to actual
20.035750</pre>
=={{header|Go}}==
Line 1,139 ⟶ 1,443:
<pre>
20.035577718385575
</pre>
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
{{works with|jq}}
'''Also works with gojq, the Go implementation of jq, and with fq'''
<syntaxhighlight lang=jq>
# output: an array
def legendreCoef($N):
{lcoef: (reduce range(0;$N+1) as $i (null; .[$i] = [range(0;$N + 1)| 0]))}
| .lcoef[0][0] = 1
| .lcoef[1][1] = 1
| reduce range(2; $N+1) as $n (.;
.lcoef[$n][0] = -($n-1) * .lcoef[$n -2][0] / $n
| reduce range (1; $n+1) as $i (.;
.lcoef[$n][$i] = ((2*$n - 1) * .lcoef[$n-1][$i-1] - ($n - 1) * .lcoef[$n-2][$i]) / $n ) )
| .lcoef ;
# input: lcoef
# output: the value
def legendreEval($n; $x):
. as $lcoef
| reduce range($n; 0 ;-1) as $i ( $lcoef[$n][$n] ; . * $x + $lcoef[$n][$i-1] ) ;
# input: lcoef
def legendreDiff($n; $x):
$n * ($x * legendreEval($n; $x) - legendreEval($n-1; $x)) / ($x*$x - 1) ;
# input: lcoef
# output: {lroots, weight}
def legendreRoots($N):
def pi: 1|atan * 4;
. as $lcoef
| { x: 0, x1: null}
| reduce range(1; 1+$N) as $i (.;
.x = ((pi * ($i - 0.25) / ($N + 0.5)) | cos )
| until (.x == .x1;
.x1 = .x
| .x as $x
| .x = .x - ($lcoef | (legendreEval($N; $x) / legendreDiff($N; $x) )) )
| .lroots[$i-1] = .x
| .x as $x
| .x1 = ($lcoef|legendreDiff($N; $x))
| .weight[$i-1] = 2 / ((1 - .x*.x) * .x1 * .x1) ) ;
# Input: {lroots, weight}
def legendreIntegrate(f; $a; $b; $N):
.lroots as $lroots
| .weight as $weight
| (($b - $a) / 2) as $c1
| (($b + $a) / 2) as $c2
| reduce range(0;$N) as $i (0; . + $weight[$i] * (($c1* $lroots[$i] + $c2)|f) )
| $c1 * .;
def task($N):
def actual: 3|exp - ((-3)|exp);
legendreCoef($N)
| legendreRoots($N)
| "Roots: ",
.lroots,
"\nWeight:",
.weight,
"\nIntegrating exp(x) over [-3, 3]: \(legendreIntegrate(exp; -3; 3; N))",
"compared to actual: \(actual)" ;
task(5)
</syntaxhighlight>
'''Invocation:'''
<pre>
jq -ncr -f gauss-legendre-quadrature.jq
</pre>
{{output}}
<pre>
Roots:
[0.906179845938664,0.5384693101056831,0,-0.5384693101056831,-0.906179845938664]
Weight:
[0.23692688505618922,0.4786286704993667,0.5688888888888889,0.4786286704993667,0.23692688505618922]
Integrating exp(x) over [-3, 3]: 20.035577718385575
compared to actual: 20.035749854819805
</pre>
Line 3,000 ⟶ 3,387:
{{trans|C}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
var N = 5
| |||