Anonymous user
Numerical integration/Gauss-Legendre Quadrature: Difference between revisions
Numerical integration/Gauss-Legendre Quadrature (view source)
Revision as of 04:42, 24 August 2018
, 5 years ago→version 3: added REXX version 3.
m (→version 2: changed wording in the REXX section header.) |
m (→version 3: added REXX version 3.) |
||
Line 2,005:
Using 82 digit precision, the N-point Gauss─Legendre quadrature (GLQ) had an accuracy of 74 digits.
</pre>
===version 3, more precision===
This REXX version is almost an exact copy of REXX version 2, but with more decimal digits of '''pi'''.
It is about twice as slow as version 2, due to the increased number of decimal digits (precision).
<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
!.= .; b= 3; a= -b; bma= b - a; bmaH= bma / 2; tiny= '1e-'digs
trueV= exp(b)-exp(a); bpa= b + a; bpaH= bpa / 2
hdr= 'iterate value (with ' digs " decimal digits being used)"
say ' step ' center(hdr, 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
/*█*/ 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 target value (Z)*/
dif= z - trueV; z= format(z, 3, digs - 2) /* " difference. */
Ndif= translate( format(dif, 3, 4, 3, 0), 'e', "E")
if #\==1 then say center(#, 6) z' ' Ndif /*no display if not computed*/
end /*#*/
say sep; xdif= compare( strip(z), trueV); say right("↑", 6 + 1 + xdif)
say left('', 6 + 1) trueV " {exact value}"; say
say 'Using ' digs " digit precision, the" ,
'N-point Gauss─Legendre quadrature (GLQ) had an accuracy of ' xdif-2 " digits."
exit /*stick a fork in it, we're all done. */
/*───────────────────────────────────────────────────────────────────────────────────────────*/
e: return 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759,
|| 4571382178525166427427466391932003059921817413596629043572900334295260595630738
/*───────────────────────────────────────────────────────────────────────────────────────────*/
pi: return 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899,
|| 862803482534211706798214808651328230664709384460955058223172535940812848111745
/*───────────────────────────────────────────────────────────────────────────────────────────*/
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= x % 1; 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
return z * e()**ix</lang>
{{out|output|text= when using the default inputs:}}
(Shown at about two-thirds size.)
<pre style="font-size:67%">
step iterate value (with 159 decimal digits being used) difference
────── ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────
2 17.4874646410555689643606840462449458421154284179349135091487247059537916662378882444064336021640614626063744948781912964250403870127054497392082425535068464109 -2.5483
3 19.8536919968055821921309108927158495960774667319753888929050027075848592516449832906645902758379575999249091274157148988582792112906526877518087112700785494497 -1.8206e-001
4 20.0286883952907008527738054439857661647073363250481518077257887668521514648379218096268747927750038360903142778646220077613647092768733641727539206268833693589 -7.0615e-003
5 20.0355777183855621539285357252750939315016272074471283081673242529514166130221254213250349496939691709537643294259047823350162410908440808868981982394287542091 -1.7214e-004
6 20.0357469750923438830654575585499253741529947892197512571761670590022501037527117346339483928363770582109285164930728028479549289382406446621705905363209981933 -2.8797e-006
7 20.0357498197266007755718729372891903369400657532378489130759167634362318526784010016150667027038415189719144094529764766032097831604495667799067330556673881546 -3.5093e-008
8 20.0357498544945172882260918041683132616236752579944055100693304551390338045262089091194019302017562870527315644307417688383478919210145963055448428522264642591 -3.2529e-010
9 20.0357498548174338368864419454858704839263168086955797931292590585320198342940085570553927472311015418220675609961921140415760514983040167737226050690228927443 -2.3700e-012
10 20.0357498548197898711175766908543458234008349625446568080936795730938134205900980645938318794902592556558231569959762420203929344018773329199723457149763574343 -1.3927e-014
11 20.0357498548198037305529147159697031241993516306485175808291929207610544866584568009626862857221858328844106864371425322111609007302709732793823163103980149653 -6.7396e-017
12 20.0357498548198037976759531014454017742327138984429607438017578771715767588391691509175808718708593063121709896967107496243434245185896147055314894150234262075 -2.7323e-019
13 20.0357498548198037979482458119092690701862659228785307035583081473361900008835808932495328864420024278695427964698380448330606714160259282675390182203803537594 -9.4143e-022
14 20.0357498548198037979491844483599375945130148356706886332919441446027039132743905494286471338717783707421873433644754993992655580745072286831502363474798175265 -2.7906e-024
15 20.0357498548198037979491872317401917248452734118643091749897281356338832738714150881537113815780435230011480697467170623887897830301712412973655748924184138648 -7.1915e-027
16 20.0357498548198037979491872389153958789316129464894982848020715833786709121310547889685984881568546203564135185474792767674806869872650180714616455691318778648 -1.6260e-029
17 20.0357498548198037979491872389316236038179252557440453906282250905385221873347716826354198555233437240574026019817833907372014036252533047705435353247648510898 -3.2517e-032
18 20.0357498548198037979491872389316560624360571301484111974244019477736095885421361807599231231543821951618639462965984321643251022835234451110049047608124949533 -5.7920e-035
19 20.0357498548198037979491872389316561202637283172074241556158972833578634894365092635000776399956033063018069653085902399896542171129596405210008317497301936898 -9.2480e-038
20 20.0357498548198037979491872389316561203560751340857503751994442223163866912408434007886096643419528065940077022083150476496426837665378721283432879108630468864 -1.3311e-040
21 20.0357498548198037979491872389316561203562080727616463861143647576984994047530870779393715057751591887673397688454357985082021265151278191050057935329724907161 -1.7360e-043
22 20.0357498548198037979491872389316561203562082461596244537077863602238433892612703628843743785373313737563806457244053157873973239947461987202443878362979981218 -2.0610e-046
23 20.0357498548198037979491872389316561203562082463655032534484950691669880046406047078766996078695370527223056578914332723730363863326194707715142045831095820995 -2.2368e-049
24 20.0357498548198037979491872389316561203562082463657267060509015976314523758814742624773428457390528961843568960502876896215809857825164102337905868347725395891 -2.2276e-052
25 20.0357498548198037979491872389316561203562082463657269286070017882824923688080311511389836619043005851350331110867389220628954338053656628671036072512306223102 -2.0430e-055
26 20.0357498548198037979491872389316561203562082463657269288111333795426189423729667519158562143832977811003145168351321839626313132075697513253761673496828204601 -1.7312e-058
27 20.0357498548198037979491872389316561203562082463657269288113063661454822050198926197665008333893008724687497228278730367375441075263700413282548634210907331356 -1.3595e-061
28 20.0357498548198037979491872389316561203562082463657269288113065019935786483820352375621786828318969009163053743757325024448325026804644277866300802833611297358 -9.9207e-065
29 20.0357498548198037979491872389316561203562082463657269288113065020927177593233999249852447888627901300469719564790181325442944469692690797774430312247159512798 -6.7451e-068
30 20.0357498548198037979491872389316561203562082463657269288113065020927851675301934062025341716601075750412806887227020916063849030412480955063639628314338766826 -4.2832e-071
31 20.0357498548198037979491872389316561203562082463657269288113065020927852103363148863217394106431702791915956948972366384835732103508918001327415359847661271744 -2.5459e-074
32 20.0357498548198037979491872389316561203562082463657269288113065020927852103617599854934274435013875248206413049448382025586066461615726348079942124364780837615 -1.4196e-077
33 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741736109635347323907131494641377410353985987829217992622815976248321170584831199 -7.4395e-081
34 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810467715704772209566910717933633388969835872983190631663850670877759345234946 -3.6713e-084
35 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504412069378446854036859408497315019337333762510854198446941961793973563786 -1.7091e-087
36 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429152646383719980280460795167918691617029439367737607466797014691070736 -7.5175e-091
37 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160160714391043273984198489693834991216803247954607301723271542659342 -3.1292e-094
38 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163842341789925349746540298990681930753381942866562579949258588756 -1.2345e-097
39 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843575809851614709658383098559963930599249691243554488598937473 -4.6221e-101
40 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576271898405984614568086424291240202255560215705599799392784 -1.6447e-104
41 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062816325200556739918251890227721352129415324255316078 -5.5685e-108
42 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062871992591098295378332977741979460337064627095382229 -1.7962e-111
43 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010547683152388008632372342584043989711962229120 -5.5262e-115
44 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553207686109280576604524821212627658539678173 -1.6234e-118
45 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309007883789778553235095713392096309190 -4.5581e-122
46 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309463568475856093621233977451367882545 -1.2245e-125
47 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309463690894669420459178304255584843389 -3.1504e-129
48 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309463690926165554457658247246849536195 -7.7695e-133
49 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309463690926173322092766217806542725001 -1.8383e-136
50 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309463690926173323930695861756195366068 -3.9547e-140
51 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309463690926173323931067751563248835441 -2.3582e-141
52 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309463690926173323931067290550361259682 -2.4043e-141
53 20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309463690926173323931250155044682676454 1.5882e-140
────── ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── ─────────────
↑
20.0357498548198037979491872389316561203562082463657269288113065020927852103617741810504429160163843576272062872010553209309463690926173323931091333618506603714 {exact value}
Using 159 digit precision, the N-point Gauss─Legendre quadrature (GLQ) had an accuracy of 141 digits.
</pre>
| |||