Anonymous user
Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Numerical integration/Gauss-Legendre Quadrature (view source)
Revision as of 23:52, 31 October 2013
, 10 years ago→{{header|REXX}}: added version 2. -- ~~~~
Walterpachl (talk | contribs) (ooRexx added (translated from REXX)) |
(→{{header|REXX}}: added version 2. -- ~~~~) |
||
Line 1,085:
=={{header|REXX}}==
===version 1===
<lang rexx>/*---------------------------------------------------------------------
* 31.10.2013 Walter Pachl Translation from PL/I
Line 1,219 ⟶ 1,220:
20.0357498548198037979491872389316561203562082463657269288113 (exact)</pre>
===version 2===
This REXX version is an optimized version (of version 1) which uses:
:::* a faster '''cos''' subroutine (which uses radian normalization)
:::* a faster '''exp''' subroutine
:::* simple variables instead of stemmed arrays
:::* static variables instead of commonly used expressions
:::* multiplication instead of division
:::* a generic approach for setting the numeric DIGITS
:::* an increase (by +1) of the number of iterations
:::* 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.*/
digs=40; d2=digs*2
numeric digits digs; tiny='1e-'d2; pi=pi(); a=-3; b=3; bma=b-a; bpa=b+a
numeric digits d2; true=exp(b)-exp(a)
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
p1z=2; p1.1=1; p1.2=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
end /*L*/
do j=1 for p1z+1;T.j=((k+km)*T.j-km*TT.j)/k; end /*j*/
p0z=p1z; do j=1 for p0z; p0.j=p1.j ; end /*j*/
p1z=p1z+1; do j=1 for p1z; p1.j= T.j ; end /*j*/
end /*k*/
do !=1 for step
x=cos(pi*(!-.25)/(step+.5))
do times%2 until abs(dx)<tiny
f=p1.1; df=0
do k=2 to p1z
df=f+x*df
f=p1.k+x*f
end /*k*/
dx=f/df; x=x-dx
end /*times%2 until···*/
r.1.!=x
r.2.!=2/((1-x*x)*df*df)
end /*!*/
sum=0
do m=1 for step
sum=sum + r.2.m * exp(bpa*.5 + r.1.m*bma*.5)
end /*m*/
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──────────────────────*/
cos: procedure; arg x; x=r2r(x); a=abs(x); numeric fuzz min(9,digits()-9)
if a=pi() then return -1; if a=pi()/2 | a=2*pi() then return 0
if a=pi()/3 then return .5; if a=2*pi()/3 then return -.5
return .sinCos(1,1,-1)
.sinCos: parse arg z 1 p,_,i; x=x*x
do k=2 by 2; _=-_*x/(k*(k+i));z=z+_;if z=p then leave;p=z;end; return z
/*──────────────────────────────────E subroutine────────────────────────*/
e: return 2.7182818284590452353602874713526624977572470936999595749669676277240766303535
/*──────────────────────────────────EXP subroutine──────────────────────*/
exp: procedure; arg x; ix=x%1; if abs(x-ix)>.5 then ix=ix+sign(x); x=x-ix
z=1; _=1; w=z; do j=1; _=_*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───────────────────────*/
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862
/*──────────────────────────────────R2R subroutine──────────────────────*/
r2r: return arg(1)//(2*pi()) /*normalize radians►1 unit circle*/</lang>
'''output'''
<pre style="overflow:scroll">
step iterative value difference
──── ─────────────────────────────────────────── ────────────
1 6.0000000000000000000000000000000000000000 -1.4036E+1
2 17.4874646410555689643606840462449458421154 -2.5483
3 19.8536919968055821921309108927158495960775 -1.8206E-1
4 20.0286883952907008527738054439857661647073 -7.0615E-3
5 20.0355777183855621539285357252750939315016 -1.7214E-4
6 20.0357469750923438830654575585499253741530 -2.8797E-6
7 20.0357498197266007755718729372891903369401 -3.5093E-8
8 20.0357498544945172882260918041683132616237 -3.2529E-10
9 20.0357498548174338368864419454858704839263 -2.3700E-12
10 20.0357498548197898711175766908543458234008 -1.3927E-14
11 20.0357498548198037305529147159697031241994 -6.7396E-17
12 20.0357498548198037976759531014454017742327 -2.7323E-19
13 20.0357498548198037979482458119092690701863 -9.4143E-22
14 20.0357498548198037979491844483599375945130 -2.7906E-24
15 20.0357498548198037979491872317401917248453 -7.1915E-27
16 20.0357498548198037979491872389153958789316 -1.6260E-29
17 20.0357498548198037979491872389316236038179 -3.2517E-32
18 20.0357498548198037979491872389316560624361 -5.7920E-35
19 20.0357498548198037979491872389316561202637 -9.2480E-38
20 20.0357498548198037979491872389316561203561 -1.3313E-40
21 20.0357498548198037979491872389316561203562 2.0044E-43
──── ─────────────────────────────────────────── ────────────
20.035749854819803797949187238931656120356208246365726928811306502092785210360717 {exact}
</pre>
=={{header|Tcl}}==
| |||