Numerical integration/Gauss-Legendre Quadrature: Difference between revisions

m
→‎version 2: added/changed comments and whitespace, changed indentations, increased precision of pi and e.
(→‎{{header|Java}}: added Java)
m (→‎version 2: added/changed comments and whitespace, changed indentations, increased precision of pi and e.)
Line 1,628:
<br>exists nested &nbsp; '''do''' &nbsp; loops with different (grouped) indentations, &nbsp; 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 pgmprogram does numerical integration using an N-point Gauss─Legendre Quadraturequadrature (GLQ)rule. */
pi=pi(); digs=length(pi)-1; numeric digits digs; reps=digs%2
!.=.; b=3; a=-b; bma=b-a; bmaH=bma/2; tiny='1E1e-' || digs
trueV=exp(b)-exp(a); bpa=b+a; bpaH=bpa/2
say ' step ' center("iterative value",digs+3) ' difference' /*show hdr*/
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
if #\==1 then say center(#, 6) z' ' Ndif /*don't display if not computed.*/
 
/*█*/ do k=2 to #; km=k-1; do y=1 for p1z; T.y=p1.y; end /*y*/
/*█*/ T.y=0; TT.=0; do L=1 for p0z; _=L+2; TT._=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 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*/
 
/*▓*/ do !=1 for #; x=cos( pi * (! - .25) / ## )
/*▓*/
/*▓*/ /*▓*/ /*░*/ do reps until abs(dx) <= tiny
/*▓*/ /*▓*/ /*░*/ f=p1.1; df=0; do u=2 to p1z
/*▓*/ /*▓*/ /*░*/ df=f + x*df
/*▓*/ /*▓*/ /*░*/ f=p1.u + x*f
/*▓*/ /*▓*/ /*░*/ end /*u*/
/*▓*/ /*▓*/ /*░*/ dx=f/df; x=x-dx
/*▓*/ /*▓*/ /*░*/ end /*reps ···*/
/*▓*/ r.1.!=x
/*▓*/ r.2.!=2 / ((1 - x**2) * df**2)
/*▓*/ end /*!*/
$=0
/*▒*/ do m=1 for #; $=$+r.2.m*exp(bpaH+r.1.m*bmaH); end /*m*/
z=bmaH*$ /*calculate the target value (Z)*/
dif=z-trueV; z=format(z, 3, digs-2) /* " /* " difference. */
Ndif=translate( format(dif, 3, 4, 2, 0), 'e', "E")
end /*#*/
 
say sep; say left('', 6+1) trueV " {exact value}"
exit /*stick a fork in it, we're all done. */
/*───────────────────────────────────────────────────────────────────────────────────────────*/
/*─────────────────────────────────────────────────────────────────────────────────────*/
e: return 2.718281828459045235360287471352662497757247093699959574966967627724076630353547595
e: return 2.7182818284590452353602874713526624977572470936999595749669676277240766303535
pi: return 3.141592653589793238462643383279502884197169399375105820974944592307816406286286209
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862
/*───────────────────────────────────────────────────────────────────────────────────────────*/
/*──────────────────────────────────────────────────────────────────────────────────*/
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=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
if z\==0 then z=z*e()**ix; return z</lang>
'''output'''
<pre>
step iterative value difference
/*────── ─────────────────────────────────────────────────────────────────────────────────────*/ ─────────────
────── ───────────────────────────────────────────────────────────────────────────────── ─────────────
2 6.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 -1.4036e+01
3 17.487464641055568964360684046244945842115428417934913509148724705953791666237348746464105556896436068404624494584211542841793491350914872470595379166623788825 -2.5483
4 19.853691996805582192130910892715849596077466731975388892905002707584859251644985369199680558219213091089271584959607746673197538889290500270758485925164498330 -1.8206e-01
5 20.028688395290700852773805443985766164707336325048151807725788766852151464837102868839529070085277380544398576616470733632504815180772578876685215146483792186 -7.0615e-03
6 20.035577718385562153928535725275093931501627207447128308167324252951416613021203557771838556215392853572527509393150162720744712830816732425295141661302212542 -1.7214e-04
7 20.035746975092343883065457558549925374152994789219751257176167059002250103751203574697509234388306545755854992537415299478921975125717616705900225010375271175 -2.8797e-06
8 20.035749819726600775571872937289190336940065753237848913075916763436231852678503574981972660077557187293728919033694006575323784891307591676343623185267840087 -3.5093e-08
9 20.035749854494517288226091804168313261623675257994405510069330455139033804525303574985449451728822609180416831326162367525799440551006933045513903380452620872 -3.2529e-10
10 20.035749854817433836886441945485870483926316808695579793129259058532019834301103574985481743383688644194548587048392631680869557979312925905853201983429400861 -2.3700e-12
11 20.035749854819789871117576690854345823400834962544656808093679573093813420607203574985481978987111757669085434582340083496254465680809367957309381342059009668 -1.3927e-14
12 20.035749854819803730552914715969703124199351630648517580829192920761054486659503574985481980373055291471596970312419935163064851758082919292076105448665845694 -6.7396e-17
13 20.035749854819803797675953101445401774232713898442960743801757877171576758880103574985481980379767595310144540177423271389844296074380175787717157675883917151 -2.7323e-19
14 20.035749854819803797948245811909269070186265922878530703558308147336190000814203574985481980379794824581190926907018626592287853070355830814733619000088357912 -9.4143e-22
15 20.035749854819803797949184448359937594513014835670688633291944144602703913361703574985481980379794918444835993759451301483567068863329194414460270391327442654 -2.7906e-24
16 20.035749854819803797949187231740191724845273411864309174989728135633883273704403574985481980379794918723174019172484527341186430917498972813563388327387142320 -7.1915e-27
17 20.035749854819803797949187238915395878931612946489498284802071583378670911519203574985481980379794918723891539587893161294648949828480207158337867091213105210 -1.6260e-29
18 20.035749854819803797949187238931623603817925255744045390628225090538522189936203574985481980379794918723893162360381792525574404539062822509053852218733547782 -3.2517e-32
19 20.035749854819803797949187238931656062436057130148411197424401947773609581645103574985481980379794918723893165606243605713014841119742440194777360958854209572 -5.7920e-35
20 20.035749854819803797949187238931656120263728317207424155615897283357863497573203574985481980379794918723893165612026372831720742415561589728335786348943623570 -9.2480e-38
21 20.035749854819803797949187238931656120356075134085750375199444222316386648510503574985481980379794918723893165612035607513408575037519944422231638669124167990 -1.3311e-40
22 20.035749854819803797949187238931656120356208072761646386114364757698499450516903574985481980379794918723893165612035620807276164638611436475769849940475037458 -1.7360e-43
23 20.035749854819803797949187238931656120356208246159624453707786360223843434844903574985481980379794918723893165612035620824615962445370778636022384338924992003 -2.0610e-46
24 20.035749854819803797949187238931656120356208246365503253448495069166988318695103574985481980379794918723893165612035620824636550325344849506916698800464997617 -2.2368e-49
25 20.035749854819803797949187238931656120356208246365726706050901597631451612071103574985481980379794918723893165612035620824636572670605090159763145237587025264 -2.2276e-52
26 20.035749854819803797949187238931656120356208246365726928607001788282492600797203574985481980379794918723893165612035620824636572692860700178828249236875179273 -2.0430e-55
27 20.035749854819803797949187238931656120356208246365726928811133379542619878657603574985481980379794918723893165612035620824636572692881113337954261894220969394 -1.7312e-58
28 20.035749854819803797949187238931656120356208246365726928811306366145483006655503574985481980379794918723893165612035620824636572692881130636614548220525870297 -1.3595e-61
29 20.03574985481980379794918723893165612035620824636572692881130650199357864896908624 -9.9207e-65
29 20.0357498548198037979491872389316561203562082463657269288113065019935767640428 -9.9208e-65
30 20.03574985481980379794918723893165612035620824636572692881130650209271775421848621 -6.7456e-68
30 20.0357498548198037979491872389316561203562082463657269288113065020927177501970 -6.7460e-68
31 20.03574985481980379794918723893165612035620824636572692881130650209278516823348154 -4.2128e-71
31 20.0357498548198037979491872389316561203562082463657269288113065020927482013630 -3.7009e-68
32 20.03574985481980379794918723893165612035620824636572692881130650209278518859457416 -2.1767e-71
32 20.0357498548198037979491872389316561203562082463657269288113065020927457504484 -3.9460e-68
────── ───────────────────────────────────────────────────────────────────────────────────── ─────────────
33 20.0357498548198037979491872389316561203562082463657269288113065020927210120532 -6.4198e-68
20.03574985481980379794918723893165612035620824636572692881130650209278521036177419 {exact value}
34 20.0357498548198037979491872389316561203562082463657269288113065020927324760748 -5.2734e-68
────── ───────────────────────────────────────────────────────────────────────────────── ─────────────
20.0357498548198037979491872389316561203562082463657269288113065020927852103607 {exact value}
</pre>