Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

Content added Content deleted
(Added Wren)
m (→‎version 2: added/changed whitespace and comments, optimized a DO loop.)
Line 2,311: Line 2,311:
<br>where there's practically no space on the right side of the REXX source statements. &nbsp; It presents a good
<br>where there's practically no space on the right side of the REXX source statements. &nbsp; It presents a good
<br>visual indication of what's what, &nbsp; but it's the dickens to pay when updating the source code.
<br>visual indication of what's what, &nbsp; but it's the dickens to pay when updating the source code.
<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)-1; 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
Line 2,319: Line 2,319:
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, 2, 0), 'e', "E")
Ndif= translate( format(dif, 3, 4, 2, 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,349: Line 2,349:
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. */
/*───────────────────────────────────────────────────────────────────────────────────────────*/
/*───────────────────────────────────────────────────────────────────────────────────────────*/
e: return 2.718281828459045235360287471352662497757247093699959574966967627724076630353547595
e: return 2.718281828459045235360287471352662497757247093699959574966967627724076630353547595
Retrieved from "https://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature"