Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

Content added Content deleted
m (→‎version 2: used more idiomatic code for computing the number of decimal digits of pi.)
m (→‎version 3, more precision: added more decimal digits for pi and e.)
Line 2,406: Line 2,406:


It is about twice as slow as version 2,   due to the doubling of the number of decimal digits   (precision).
It is about twice as slow as version 2,   due to the doubling of the number of decimal digits   (precision).
<lang rexx>/*REXX program does numerical integration using an N-point Gauss─Legendre quadrature rule. */
<lang rexx>/*REXX program does numerical integration using an N─point Gauss─Legendre quadrature rule. */
pi= pi(); digs= length(pi)-1; numeric digits digs; reps= digs % 2
pi= pi(); digs= length(pi) - length(.); numeric digits digs; reps= digs % 2
!.= .; b= 3; a= -b; bma= b - a; bmaH= bma / 2; tiny= '1e-'digs
!.= .; b= 3; a= -b; bma= b - a; bmaH= bma / 2; tiny= '1e-'digs
trueV= exp(b)-exp(a); bpa= b + a; bpaH= bpa / 2
trueV= exp(b)-exp(a); bpa= b + a; bpaH= bpa / 2
Line 2,414: Line 2,414:
sep='──────' copies("─", digs+3) '─────────────'; say sep /* " sep*/
sep='──────' copies("─", digs+3) '─────────────'; say sep /* " sep*/


do #=1 until dif>0; p0z= 1; p0.1= 1; p1z= 2; p1.1= 1; p1.2= 0; ##= # + .5; r.= 0
do #=1 until dif>0; p0z= 1; p0.1= 1; p1z= 2; p1.1= 1; p1.2= 0; ##= # + .5; r.= 0
/*█*/ do k=2 to #; km= k - 1; do y=1 for p1z; T.y= p1.y; end /*y*/
/*█*/ do k=2 to #; km= k - 1
/*█*/ T.y= 0; TT.= 0; do L=1 for p0z; _= L + 2; TT._= p0.L; end /*L*/
/*█*/ do y=1 for p1z; T.y= p1.y; end /*y*/
/*█*/ kkm= k + km; do j=1 for p1z +1; T.j= (kkm*T.j -km*TT.j)/k; end /*j*/
/*█*/ T.y= 0; TT.= 0; do L=1 for p0z; _= L + 2; TT._= p0.L; end /*L*/
/*█*/ p0z= p1z; do n=1 for p0z; p0.n= p1.n ; end /*n*/
/*█*/ kkm= k + km; do j=1 for p1z +1; T.j= (kkm*T.j - km*TT.j)/k; end /*j*/
/*█*/ p1z= p1z + 1; do p=1 for p1z; p1.p= T.p ; end /*p*/
/*█*/ p0z= p1z; do n=1 for p0z; p0.n= p1.n ; end /*n*/
/*█*/ p1z= p1z + 1; do p=1 for p1z; p1.p= T.p ; end /*p*/
/*█*/ end /*k*/
/*█*/ end /*k*/
/*▓*/ do !=1 for #; x= cos( pi * (! - .25) / ## )
/*▓*/ do !=1 for #; x= cos( pi * (! - .25) / ## )
/*▓*/
/*▓*/ /*░*/ do reps until abs(dx) <= tiny
/*▓*/ /*░*/ do reps until abs(dx) <= tiny
/*▓*/ /*░*/ f= p1.1; df= 0; do u=2 to p1z; df= f + x*df
/*▓*/ /*░*/ f= p1.1; df= 0; do u=2 to p1z; df= f + x*df
/*▓*/ /*░*/ f= p1.u +x*f
/*▓*/ /*░*/ f= p1.u +x*f
/*▓*/ /*░*/ end /*u*/
/*▓*/ /*░*/ end /*u*/
/*▓*/ /*░*/ dx= f / df; x= x - dx
/*▓*/ /*░*/ dx= f / df; x= x - dx
/*▓*/ /*░*/ end /*reps ···*/
/*▓*/ /*░*/ end /*reps ···*/
/*▓*/ r.1.!= x
/*▓*/ r.1.!= x
/*▓*/ r.2.!= 2 / ( (1 - x*x) * df*df)
/*▓*/ r.2.!= 2 / ( (1 - x**2) * df**2)
/*▓*/ end /*!*/
/*▓*/ end /*!*/
$= 0
$= 0
/*▒*/ do m=1 for #; $=$ + r.2.m * exp(bpaH + r.1.m*bmaH); end /*m*/
/*▒*/ do m=1 for #; $=$ + r.2.m * exp(bpaH + r.1.m*bmaH); end /*m*/
z= bmaH * $ /*calculate target value (Z)*/
z= bmaH * $ /*calculate target value (Z)*/
dif= z - trueV; z= format(z, 3, digs - 2) /* " difference. */
dif= z - trueV; z= format(z, 3, digs - 2) /* " difference. */
Ndif= translate( format(dif, 3, 4, 3, 0), 'e', "E")
Ndif= translate( format(dif, 3, 4, 3, 0), 'e', "E")
if #\==1 then say center(#, 6) z' ' Ndif /*no display if not computed*/
if #\==1 then say center(#, 6) z' ' Ndif /*no display if not computed*/
Line 2,444: Line 2,444:
say 'Using ' digs " digit precision, the" ,
say 'Using ' digs " digit precision, the" ,
'N-point Gauss─Legendre quadrature (GLQ) had an accuracy of ' xdif-2 " digits."
'N-point Gauss─Legendre quadrature (GLQ) had an accuracy of ' xdif-2 " digits."
exit /*stick a fork in it, we're all done. */
exit 0 /*stick a fork in it, we're all done. */
/*───────────────────────────────────────────────────────────────────────────────────────────*/
cos: procedure expose !.; parse arg x; if !.x\==. then return !.x; _= 1; z=1; y= x*x
do k=2 by 2 until p==z; p=z; _= -_*y/(k*(k-1)); z=z+_; end; !.x=z; return z
/*───────────────────────────────────────────────────────────────────────────────────────────*/
exp: procedure; parse arg x; ix= x % 1; if abs(x-ix)>.5 then ix= ix + sign(x); x= x-ix; z= 1
_=1; do j=1 until p==z; p=z; _= _*x/j; z= z+_; end; return z * e()**ix
/*───────────────────────────────────────────────────────────────────────────────────────────*/
/*───────────────────────────────────────────────────────────────────────────────────────────*/

e: return 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759,
e: return 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759,
||457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794
|| 4571382178525166427427466391932003059921817413596629043572900334295260595630738
/*───────────────────────────────────────────────────────────────────────────────────────────*/
/*───────────────────────────────────────────────────────────────────────────────────────────*/
pi: return 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899,
pi: return 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899,
||862803482534211706798214808651328230664709384460955058223172535940812848111745028410270194
|| 8628034825342117067982148086513282306647093844609550582231725359408128481117450
/*───────────────────────────────────────────────────────────────────────────────────────────*/
/*───────────────────────────────────────────────────────────────────────────────────────────*/
cos: procedure expose !.; parse arg x; if !.x\==. then return !.x; _= 1; z=1; y= x*x
cos: procedure expose !.; parse arg x; if !.x\==. then return !.x; _= 1; z=1; y= x*x
Retrieved from "https://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature"