Rosetta Code

Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

Page Discussion

Numerical integration/Gauss-Legendre Quadrature (view source)

Revision as of 21:46, 1 November 2013

117 bytes added ,  10 years ago
m
→‎version 2: re-instated the latest version for REXX version 2. -- ~~~~
m (→‎version 1: corrected time 'E' to time 'R' (reset))
m (→‎version 2: re-instated the latest version for REXX version 2. -- ~~~~)
Line 1,235:
:::*   about 1/6 the computing time
Note that the function values for   '''pi'''   and   '''e'''   should have more precision than the number of digits specified.
<lang rexx>/*REXX pgm does numerical integration using Gauss-Legendre Quadrature. */
call time 'Reset'
digs=40; d2=digs*2
numericdigs=40; digits d2=digs*2; tiny = '1e1E-'d2; pi=pi(); a=-3; b=3; bma=b-a;/*equivalent bpa=b+ato 1÷10**80 */
numeric digits d2digs; true pi=exppi(b); a=-3; b=3; bma =b-exp(a); bpa =b+a
numeric digits digs+10d2; true = exp(b)-exp(a); bmaH=bma/2; bpaH=bpa/2
numeric digits digs+10; times = digs%2 + 1
say 'step' center("iterative value",digs+3) ' difference ' /*show hdr*/
sep='────' copies('─' ,digs+3) '────────────'; say sep
 
do step=1 for times; r.=0; p0z=1; p0.1=1; step_=step+.5
p1z=2; p1.1=1; p1.2=0; r.=0
/*█*/ do k=2 to step; km=k-1
/*█*/ do L=1 for p1z; T.L=p1.L
/*█*/ end /*L*/
/*█*/ T.L=0; TT.=0
/*█*/ do L=1 for p0z; L2=L+2; TT.L2=p0.L
/*█*/ TT.L2=p0.end /*L*/
/*█*/
end /*L*/▼
/*█*/ p1z=p1z+1; do j=1 for p1z+1; p1T.j= ((k+km)*T.j -km*TT.j)/k; end /*j*/▼
 
/*█*/ p0z=p1z; do j=1 for p1z+1p0z;T p0.j=((k+km)*Tp1.j-km*TT.j)/k ; end /*j*/
/*█*/ p0zp1z=p1z+1; do j=1 for p0zp1z; p0p1.j=p1 T.j ; end /*j*/
/*█*/ end /*k*/▼
▲ p1z=p1z+1; do j=1 for p1z; p1.j= T.j ; end /*j*/
▲ do !=1 for end /*L*/step
▲ end /*k*/
r.2.! x=2/cos(pi*(1!-x*x.25)*df*df/step_)▼
 
do !=1times%2 foruntil stepabs(dx)<tiny
x f=cos(pi*(!-p1.25)/(step+.5))1; df=0
do times%2 until abs(dx)<tiny do k=2 to p1z
f=p1.1; df=0f+x*df
do f=p1.k=2 to p1z+x*f
df=f+x end /*dfk*/
dx=f/df; fx=p1.k+x*f-dx
end /*ktimes%2 until···*/
dxr.1.!=f/df; x=x-dx
end r.2.!=2/((1-x*times%2 until···x)*/df*df)
r.1. end /*!=x*/
▲ r.2.!=2/((1-x*x)*df*df)
end /*!*/
sum=0
do m=1 for step; sum=sum+r.2.m*exp(bpaH+r.1.m*bmaH)
sum=sum end + r.2.m /* exp(bpa*.5 + r.1.m*bma*.5)/
z=bmabmaH*sum*.5▼
end /*m*/
say right(step,4) format(z,2,digs) translate(format(z-true,2,4,,0),'e',"E")▼
▲ z=bma*sum*.5
▲ say right(step,4) format(z,2,digs) format(z-true,2,4,,0)
end /*step*/
 
say sep; say left('',4) true " {exact}"
say '... and took' format(time('e'),,2) "seconds."
exit /*stick a fork in it, we're done.*/
/*──────────────────────────────────COS subroutine───────────────────────────────subroutine──────────────────────*/
cos: procedure; x=r2r(arg(1)); numeric fuzz min(9,digits()-9); _=1; z=1; p=1; x=x*x
y=x*x; do k=2 by 2 until p=z; p=z; _=-_*xy/(k*(k-1)); z=z+_; if z=p thenend; return z; p=z; end
/*──────────────────────────────────E subroutine──────────────────────────────────────*/
e: return 2.7182818284590452353602874713526624977572470936999595749669676277240766303535
/*──────────────────────────────────EXP subroutine────────────────────────────────────subroutine──────────────────────*/
exp: procedure; parse arg x; ix=x%1; if abs(x-ix)>.5 then ix=ix+sign(x); x=x-ix
x=x-ix; z=1; _=1; w=z; do j=1 until p==z; p=z; _=_*x/j; z=(z+_)/1; if z==w then leave; w=z; end
if z\==0 then z=z*e()**ix; return z
/*──────────────────────────────────PI subroutine──────────────────────────────────────*/
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/183733"