Rosetta Code

Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

Page Discussion

Numerical integration/Gauss-Legendre Quadrature (view source)

Revision as of 08:00, 25 September 2015

30 bytes removed ,  8 years ago
→‎version 2: added/changed whitespace and comments, changed output alignment and indentation, used a better method of determining the numeric (decimal digits) precision.
(→‎version 2: added/changed whitespace and comments, changed output alignment and indentation, used a better method of determining the numeric (decimal digits) precision.)
Line 1,378:
 
===version 2===
This REXX version is (an optimized version (of version 1) which  and uses:
:::*   a faster   '''cos'''   function   (with full precision)
:::*   a faster   '''exp'''   function   (with full precision)
Line 1,384:
:::*   some static variables instead of repeated expressions
:::*   calculations using full (specified) precision (''numeric digits'')
:::*   multiplication using   [···'''*.5''']   instead of division using   [···'''/2''']
:::*   a generic approach for setting the   ''numeric digits''
:::*   a better test for earlier termination (stopping) of calculations
<br><br>The use of "vertical bars" is one of the very few times to use leading comments, as there isn't that many situations where there ▼
Note that the function values for &nbsp; '''pi''' &nbsp; and &nbsp; '''e''' &nbsp; should have more precision than the number of digits specified.
<br>exists nested &nbsp; DO'''do''' &nbsp; loops with different (grouped) indentations, &nbsp; and practically no space on the right side of the statements.▼
▲<br><br>The use of vertical bars is one of the very few times to use leading comments, as there isn't that many situations where there
▲<br>exists nested &nbsp; DO &nbsp; loops with different (grouped) indentations, and practically no space on the right side of the statements.
<br>It presents a good visual indication of what's what, but it's the dickens to pay when updating the code.
<lang rexx>/*REXX pgm doesperforms numerical integration using GLQ: Gauss─Legendre Quadrature. */
parsepi=pi(); arg digs .=length(pi); ifnumeric digits digs==''; then digs=70 /*assume the DIGS default?*/times=digs%2
numeric!.=.; digits digs*2+5 b=3; a=-b; bma=b-a; bmaH=bma/2; /*usetiny='1e-' higher|| working DIGs.*/digits()
timestrueV=digs%2exp(b)-exp(a); b=3; a=-b; bma bpa=b-+a; bmaH bpaH=bmabpa/2; tiny oldD='1E-' || (digs*2)8e8
say ' step ' center("iterative value",digs+53) ' difference' /*show hdr*/▼
trueV=exp(b)-exp(a); bpa=b+a; bpaH=bpa/2; oldZ=
sep='──────' copies('─' ,digs+53) '─────────────'; say sep▼
numeric digits digs+15; pi=pi(); !.=. /*use lower working DIGITs*/
▲say ' step ' center("iterative value",digs+5) ' difference' /*show hdr*/
▲sep='──────' copies('─' ,digs+5) '─────────────'; say sep
 
do step=1; for times; p0z=1; p0.1=1; step_=step + .5; r.=0
p1z=2; p1.1=1; p1.2=0; r.=0
/*█*/ do k=2 to step; km=k-1
/*█*/ do Ly=1 for p1z; T.Ly=p1.Ly; end /*y*/
/*█*/ T.y=0; end /*L*/TT.=0
/*█*/ T. do L=01 for p0z; _=L+2; TT._=0p0.L; end /*L*/
/*█*/ do L=1 for p0z; L2=L+2; TT.L2=p0.L
/*█*/ end /*L*/
/*█*/
/*█*/ kkm=k+km; do j=1 for p1z+1; T.j=(kkm*T.j-km*TT.j)/k ; end /*j*/
/*█*/ p0z=p1z; do jn=1 for p0z; p0.jn=p1.jn ; end /*jn*/
/*█*/ p1z=p1z+1; do jp=1 for p1z; p1.jp= T.jp ; end /*jp*/
/*█*/ end /*k*/
 
/*▓*/ do !=1 for step
/*▓*/ x=cos(pi*(!-.25)/step_)
/*▓*/ do times%2 until abs(dx) < tiny
/*▓*/ f=p1.1; df=0; do k=2 to p1z
/*▓*/ do k df=2f to+ p1zx*df
/*▓*/ df=f f=p1.k + x*dff
/*▓*/ f=p1.k + x end /*k*f/
/*▓*/ dx=f/df; end /*k*/x=x-dx
/*▓*/ dx=f/df; x=x-dx end /*times ···*/
/*▓*/ end /*times%2 ···*/
/*▓*/ r.1.!=x
/*▓*/ r.2.!=2 / ((1 - x*x*2) * df * df*2)
/*▓*/ end /*!*/
sum=0
/*▒*/ do m=1 for step; sum=sum + r.2.m * exp(bpaH+r.1.m*bmaH)
/*▒*/ end /*m*/
z=bmaH*sum /*calculate the target value. (Z)*/
dif=translate(z-trueV,; 'e', "E" z=format(z,3,digs-2) /* " /*calculate the " difference. */
if abs(dif)>=oldD then leave; Ndif=translate( format(dif,3,4,2,0), 'e', "E")
z=format(z,3,digs); if z==oldZ then leave /*exhausted computations? */
say center(step,6) z' ' format(dif,3,4,2,0) Ndif /*display target and diffdifference.*/
oldD=abs(dif)
oldZ=z /*new target for the oldZ.*/
end /*step*/ /* [↑] Z is differentthe difference. */
 
say sep; say left('',6+1) format( trueV) " {exact value}"
exit /*stick a fork in it, we're all done. */
/*─────────────────────────────────────────────────────────────────────────────────────*/
/*──────────────────────────────────one-liner subroutines───────────────*/
e: return 2.7182818284590452353602874713526624977572470936999595749669676277240766303535
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862
/*──────────────────────────────────────────────────────────────────────────────────*/
/*──────────────────────────────────COS subroutine──────────────────────────────────*/
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 subroutine────────────────────────────*/
exp: procedure; parse arg x; ix=trunc(x); 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
Line 1,451 ⟶ 1,445:
'''output'''
<pre>
step iterative value difference
────── ───────────────────────────────────────────────────────────────────────────────── ─────────────
────── ─────────────────────────────────────────────────────────────────────────── ─────────────
1 6.0000000000000000000000000000000000000000000000000000000000000000000000000000 -1.4036e+01
1 6.0000000000000000000000000000000000000000000000000000000000000000000000 -1.4036E+01
2 17.48746464105556896436068404624494584211542841793491350914872470595379174874646410555689643606840462449458421154284179349135091487247059537916662373 -2.5483
3 19.8536919968055821921309108927158495960774667319753888929050027075848592516449 -1.8206e-01
3 19.8536919968055821921309108927158495960774667319753888929050027075848593 -1.8206E-01
4 20.0286883952907008527738054439857661647073363250481518077257887668521514648371 -7.0615e-03
4 20.0286883952907008527738054439857661647073363250481518077257887668521515 -7.0615E-03
5 20.0355777183855621539285357252750939315016272074471283081673242529514166130212 -1.7214e-04
5 20.0355777183855621539285357252750939315016272074471283081673242529514166 -1.7214E-04
6 20.0357469750923438830654575585499253741529947892197512571761670590022501037512 -2.8797e-06
6 20.0357469750923438830654575585499253741529947892197512571761670590022501 -2.8797E-06
7 20.0357498197266007755718729372891903369400657532378489130759167634362318526785 -3.5093e-08
7 20.0357498197266007755718729372891903369400657532378489130759167634362319 -3.5093E-08
8 20.0357498544945172882260918041683132616236752579944055100693304551390338045253 -3.2529e-10
8 20.0357498544945172882260918041683132616236752579944055100693304551390338 -3.2529E-10
9 20.0357498548174338368864419454858704839263168086955797931292590585320198343011 -2.3700e-12
9 20.0357498548174338368864419454858704839263168086955797931292590585320198 -2.3700E-12
10 20.0357498548197898711175766908543458234008349625446568080936795730938134206072 -1.3927e-14
10 20.0357498548197898711175766908543458234008349625446568080936795730938134 -1.3927E-14
11 20.0357498548198037305529147159697031241993516306485175808291929207610544866595 -6.7396e-17
11 20.0357498548198037305529147159697031241993516306485175808291929207610545 -6.7396E-17
12 20.0357498548198037976759531014454017742327138984429607438017578771715767588801 -2.7323e-19
12 20.0357498548198037976759531014454017742327138984429607438017578771715768 -2.7323E-19
13 20.0357498548198037979482458119092690701862659228785307035583081473361900008142 -9.4143e-22
13 20.0357498548198037979482458119092690701862659228785307035583081473361900 -9.4143E-22
14 20.0357498548198037979491844483599375945130148356706886332919441446027039133617 -2.7906e-24
14 20.0357498548198037979491844483599375945130148356706886332919441446027039 -2.7906E-24
15 20.0357498548198037979491872317401917248452734118643091749897281356338832737044 -7.1915e-27
15 20.0357498548198037979491872317401917248452734118643091749897281356338833 -7.1915E-27
16 20.0357498548198037979491872389153958789316129464894982848020715833786709115192 -1.6260e-29
16 20.0357498548198037979491872389153958789316129464894982848020715833786709 -1.6260E-29
17 20.0357498548198037979491872389316236038179252557440453906282250905385221899362 -3.2517e-32
17 20.0357498548198037979491872389316236038179252557440453906282250905385222 -3.2517E-32
18 20.0357498548198037979491872389316560624360571301484111974244019477736095816451 -5.7920e-35
18 20.0357498548198037979491872389316560624360571301484111974244019477736096 -5.7920E-35
19 20.0357498548198037979491872389316561202637283172074241556158972833578634975732 -9.2480e-38
19 20.0357498548198037979491872389316561202637283172074241556158972833578635 -9.2480E-38
20 20.0357498548198037979491872389316561203560751340857503751994442223163866485105 -1.3311e-40
20 20.0357498548198037979491872389316561203560751340857503751994442223163867 -1.3311E-40
21 20.0357498548198037979491872389316561203562080727616463861143647576984994505169 -1.7360e-43
21 20.0357498548198037979491872389316561203562080727616463861143647576984994 -1.7360E-43
22 20.0357498548198037979491872389316561203562082461596244537077863602238434348449 -2.0610e-46
22 20.0357498548198037979491872389316561203562082461596244537077863602238434 -2.0610E-46
23 20.0357498548198037979491872389316561203562082463655032534484950691669883186951 -2.2368e-49
23 20.0357498548198037979491872389316561203562082463655032534484950691669880 -2.2368E-49
24 20.0357498548198037979491872389316561203562082463657267060509015976314516120711 -2.2276e-52
24 20.0357498548198037979491872389316561203562082463657267060509015976314524 -2.2276E-52
25 20.0357498548198037979491872389316561203562082463657269286070017882824926007972 -2.0430e-55
25 20.0357498548198037979491872389316561203562082463657269286070017882824924 -2.0430E-55
26 20.0357498548198037979491872389316561203562082463657269288111333795426198786576 -1.7312e-58
26 20.0357498548198037979491872389316561203562082463657269288111333795426189 -1.7312E-58
27 20.0357498548198037979491872389316561203562082463657269288113063661454830066555 -1.3595e-61
27 20.0357498548198037979491872389316561203562082463657269288113063661454822 -1.3595E-61
28 20.0357498548198037979491872389316561203562082463657269288113065019935767640428 -9.9208e-65
28 20.0357498548198037979491872389316561203562082463657269288113065019935786 -9.9207E-65
29 20.0357498548198037979491872389316561203562082463657269288113065020927177501970 -6.7460e-68
29 20.0357498548198037979491872389316561203562082463657269288113065020927178 -6.7451E-68
30 20.0357498548198037979491872389316561203562082463657269288113065020927482013630 -3.7009e-68
30 20.0357498548198037979491872389316561203562082463657269288113065020927852 -4.2829E-71
────── ───────────────────────────────────────────────────────────────────────────────── ─────────────
────── ─────────────────────────────────────────────────────────────────────────── ─────────────
20.035749854819803797949187238931656120356208246365726928811306502092785210360716529000357498548198037979491872389316561203562082463657269288113065020927852103607 {exact value}
</pre>
 
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/183764"