Difference between revisions of "Integral"
Line 25:  Line 25:  
==Introduction==  ==Introduction==  
−  +  According to [https://en.wikipedia.org/wiki/Numerical_integration Wikipedia], "In analysis, '''numerical integration''' comprises a broad family of algorithms for calculating the numerical value of a definite integral".  
−  
−  
−  
−  
+  This implementation uses several different algorithms to achieve this as well as extends these algorithms to Infinite and Left and Right semiinfinite intervals.  
==Notation==  ==Notation==  
Line 73:  Line 70:  
==Variants==  ==Variants==  
−  There are several different algorithms that may be used for Numerical Integration  +  There are several different algorithms that may be used for Numerical Integration; the main (and default) one is based upon the TanhSinh Quadrature code in David H. Bailey's [https://www.davidhbailey.com/dhbsoftware/arprec2.2.18.tar.gz ARPREC] package. Other algorithms include [https://en.wikipedia.org/wiki/Gaussian_quadrature GaussLegendre] and [https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas NewtonCotes] taken from Laurent Fousse's CRQ (Correctly Rounded Quadrature) [https://github.com/lfousse/crq Numerical Integration] code written in MPFR as described in [https://tel.archivesouvertes.fr/tel00477243/document this paper]. These algorithm may be selected via the Variant operator as in 
<apll><pre>  <apll><pre>  
Line 80:  Line 77:  
6.2831853071795864769252867665590057683<_mark>9</_mark>  6.2831853071795864769252867665590057683<_mark>9</_mark>  
{1+1○⍵}∫⍠'n' ○2<_x/> ⍝ NewtonCotes  {1+1○⍵}∫⍠'n' ○2<_x/> ⍝ NewtonCotes  
−  6.  +  6.283185307179586476925286766559005<_mark>7</_mark>1385 
○2<_x/> ⍝ Exact answer  ○2<_x/> ⍝ Exact answer  
6.28318530717958647692528676655900576839  6.28318530717958647692528676655900576839  
</pre></apll>  </pre></apll>  
−  The Order of the Numerical Integration (the number of rectangles used to approximate the result) is, by default, <apll>128</apll> for the GaussLegendre algorithm, and <apll>  +  The Order of the Numerical Integration (the number of rectangles used to approximate the result) is, by default, <apll>128</apll> for the GaussLegendre algorithm, and <apll>32</apll> for NewtonCotes. That number may be changed via the Variant operator as in 
<apll><pre>  <apll><pre>  
Line 91:  Line 88:  
⎕PP←60 ⋄ ⎕FPC←512  ⎕PP←60 ⋄ ⎕FPC←512  
{÷1+⍵*2}∫⍠'g' N  {÷1+⍵*2}∫⍠'g' N  
−  1.  +  1.5120405040791739263291383891879796566219342815871082643439<_mark>7</_mark> 
−  {÷1+⍵*2}∫⍠('  +  {÷1+⍵*2}∫⍠'n' N 
−  1.  +  1.512040504079173926329138<_mark>3</_mark>4207540700696171248671126465236073 
−  {÷1+⍵*2}∫⍠('  +  {÷1+⍵*2}∫⍠('n' 64) N 
−  1.  +  1.51204050407917392632913838918797965662193428158<_mark>7</_mark>22152948792 
+  {÷1+⍵*2}∫⍠('n' 96) N  
+  1.5120405040791739263291383891879796566219342815871082643439<_mark>7</_mark>  
¯3○N ⍝ Exact answer  ¯3○N ⍝ Exact answer  
−  
</pre></apll>  </pre></apll>  
+  
+  <p>There are additional choices for algorithms such as</p>  
+  
+  <table border="1" cellpadding="5" cellspacing="0" rules="none" summary="">  
+  <tr>  
+  <th align="left">Algorithm</th>  
+  <th align="left">Syntax</th>  
+  <th align="left">Source</th>  
+  </tr>  
+  
+  <tr>  
+  <td>GaussLegendre</td>  
+  <td><apll><_sg/>⍠'G'</apll></td>  
+  <td>David H. Bailey's ARPREC ARPREC C++ code</td>  
+  </tr>  
+  
+  <tr>  
+  <td>Error Function</td>  
+  <td><apll><_sg/>⍠'E'</apll></td>  
+  <td>David H. Bailey's ARPREC ARPREC C++ code</td>  
+  </tr>  
+  
+  <tr>  
+  <td>TanhSinh</td>  
+  <td><apll><_sg/>⍠'T'</apll></td>  
+  <td>David H. Bailey's ARPREC ARPREC C++ code</td>  
+  </tr>  
+  
+  <tr>  
+  <td valign="top">TanhSinh</td>  
+  <td><apll><_sg/>⍠'t'</apll> and <apll><_sg/></apll><br/>as it is the default algorithm</td>  
+  <td valign="top">David H. Bailey's ARPREC ARPREC C++ code translated to C</td>  
+  </tr>  
+  
+  <tr>  
+  <td>GaussLegendre</td>  
+  <td><apll><_sg/>⍠'g'</apll></td>  
+  <td>Laurent Fousse's CRQ code</td>  
+  </tr>  
+  
+  <tr>  
+  <td>NewtonCotes</td>  
+  <td><apll><_sg/>⍠'n'</apll></td>  
+  <td>Laurent Fousse's CRQ code</td>  
+  </tr>  
+  
+  <tr>  
+  <td>TanhSinh</td>  
+  <td><apll><_sg/>⍠'d'</apll></td>  
+  <td>Graeme Dennes Excel Spreadsheet code</td>  
+  </tr>  
+  </table>  
+  
+  <p>When the Variant operator is applied to the default algorithm, the first number is used to replace the exponent <apll>N=5</apll> in the loop limit of <apll>12×2*N</apll>.</p>  
+  
+  <p>Not all of these additional algorithm's may be supported in future versions.</p>  
==Examples==  ==Examples==  
−  +  The formula for the [https://en.wikipedia.org/wiki/Parabola#Similarity_to_the_unit_parabola Unit Parabola] is <apll>{⍵*2}</apll> whose Integral is <apll>{1<_r/>3×⍵*3}</apll>, and so the Integral of the Unit Parabola from <apll>a</apll> to <apll>a</apll> is <apll>/{1<_r/>3×⍵*3}a,a ←→ 2<_r/>3×a*3</apll>:  
−  +  <apll><pre>  
+  ⎕FPC←128  
+  a←2 ⋄ (a){⍵*2}∫a  
+  5.33333333333333333333333333333333333334  
+  </pre></apll>  
+  
+  This may be used to calculate the [https://en.wikipedia.org/wiki/Parabola#Area_enclosed_between_a_parabola_and_a_chord Area Enclosed Between a Parabola and a Chord], the formula for which is twothirds of the area of the surrounding rectangle of the chord, the line tangent to the Vertex, and the two sides <apll>x=a</apll> and <apll>x=a</apll>.  
+  
+  Because the Unit Parabola curves upward and away from the Xaxis, the above Integral returns the area under (i.e. below and outside) the Parabola. In this case, because we're Integrating from <apll>a</apll> to <apll>a</apll>, the chord is the line <apll>y={⍵*2} a ←→ y=a*2</apll>. The top of the rectangle (chord) is of length <apll>aa ←→ 2×a</apll> and of height <apll>y ←→ a*2</apll>, whose area is <apll>(2×a)×a*2 ←→ 2×a*3</apll>. Subtracting the result of the Integral from the rectangle's area and multiplying this by twothirds yields the area enclosed between a chord and a Parabola.  
<apll><pre>  <apll><pre>  
−  +  2<_r/>3×(2×a*3)(a){⍵*2}∫a  
−  +  7.1111111111111111111111111111111111111  
−  
</pre></apll>  </pre></apll>  
Line 115:  Line 176:  
<apll><pre>  <apll><pre>  
−  {1+1○⍵}  +  {1+1○⍵}<_sg/>○2<_x/> ⋄ ○2<_x/> 
−  6.  +  6.283185307179586476925286766559005768<_mark>3</_mark>4 
6.28318530717958647692528676655900576839  6.28318530717958647692528676655900576839  
</pre></apll>  </pre></apll>  
Line 124:  Line 185:  
<apll><pre>  <apll><pre>  
{1○⍵}∫○2<_x/>  {1○⍵}∫○2<_x/>  
−  +  6.69453502284549420282962009846618284241<_E/>¯42  
</pre></apll>  </pre></apll>  
−  A [https://en.wikipedia.org/wiki/Normal_distribution Normal Distribution] is defined as <apll>nd←{(*¯0.5×⍵*2)÷√○2<_x/>}</apll>.  +  A [https://en.wikipedia.org/wiki/Normal_distribution Normal Distribution] is defined as <apll>nd←{(*¯0.5×⍵*2)÷√○2<_x/>}</apll>. Because this instance of the algorithm is scaled properly, integrating it over the entire width from <apll>¯∞</apll> to <apll>∞</apll> yields an answer of <apll>1</apll> (the area under the curve). However, to achieve this result, we need to increase the default loop limit exponent: 
<apll><pre>  <apll><pre>  
−  +  ¯∞ nd∫ ∞  
−  1.  +  1.000000000000002522695227400560226259006 
+  ¯∞ nd∫⍠10 ∞  
+  1  
</pre></apll>  </pre></apll>  
Line 138:  Line 201:  
<apll><pre>  <apll><pre>  
⍪¯1 ¯2 ¯3 nd∫¨ 1 2 3  ⍪¯1 ¯2 ¯3 nd∫¨ 1 2 3  
−  0.  +  0.682689492137085897170465091264075844958 ⍝ 68% 
−  0.  +  0.954499736103641585599434725666933125059 ⍝ 95% 
0.997300203936739810946696370464810045244 ⍝ 99.7%  0.997300203936739810946696370464810045244 ⍝ 99.7%  
</pre></apll>  </pre></apll>  
Line 147:  Line 210:  
==Numerical Differentiation==  ==Numerical Differentiation==  
−  Note that [[Derivative]], the inverse of this operator, has also been implemented.  +  Note that [[Derivative]], the inverse of this operator, has also been implemented. These two features may be combined in various ways, for example, to calculate [https://en.wikipedia.org/wiki/Arc_length#Finding_arc_lengths_by_integrating Arc Length]. 
−  +  The formula for Arc Length is simply <apll>arclen←{4○⍺⍺∂⍵}</apll> which can be used on any given function. For example, the upper half of a unit circle is represented by the function <apll>{0○⍵}</apll>. The interval between <apll>[√0.5, √0.5]</apll> is then a quarter of a circle, whose length is  
−  +  <apll><pre>  
+  ⎕FPC←512 ⋄ ⎕PP←60  
+  (√0.5<_x/>){0○⍵}arclen∫√0.5<_x/>  
+  1.5707963267948966192313216916397514420985846996875529104874<_mark>7</_mark>  
+  ○0.5<_x/> ⍝ Exact length  
+  1.57079632679489661923132169163975144209858469968755291048747</pre></apll> 
Revision as of 17:09, 27 March 2020


L is an optional Real numeric singleton which represents the lower bound of the definite Integral. If it is omitted, 0 is used.  
R is a Real numeric singleton which represents the upper bound of the definite Integral.  
f is an arbitrary monadic function whose argument and result are both Real numeric singletons. 
Introduction
According to Wikipedia, "In analysis, numerical integration comprises a broad family of algorithms for calculating the numerical value of a definite integral".
This implementation uses several different algorithms to achieve this as well as extends these algorithms to Infinite and Left and Right semiinfinite intervals.
Notation
The symbol chosen for this operator is the Integral Sign (∫), entered from the keyboard as Alt’S’ or AltShift's' (U+222B), and used in mathematics for Integration.
Applications
There are many, many applications of Numerical Integration. Here are but a few taken from Whitman College's Mathematics Department's Calculus Online course:
 Area between curves
 Distance, Velocity, Acceleration
 Volume
 Average value of a function
 Work
 Center of Mass
 Kinetic energy; improper integrals
 Probability
 Arc Length
 Surface Area
and more taken from the website Interactive Mathematics:
 Applications of the Indefinite Integral
 Area Under a Curve by Integration
 Area Between 2 Curves using Integration
 Volume of Solid of Revolution by Integration
 Shell Method: Volume of Solid of Revolution
 Centroid of an Area by Integration
 Moments of Inertia by Integration
 Work by a Variable Force using Integration</a>
 Electric Charges by Integration
 Average Value of a Function by Integration
 Force Due to Liquid Pressure by Integration</a>
 Arc Length of a Curve using Integration</a>
 Arc Length of Curve: Parametric, Polar Coordinates</a>
Variants
There are several different algorithms that may be used for Numerical Integration; the main (and default) one is based upon the TanhSinh Quadrature code in David H. Bailey's ARPREC package. Other algorithms include GaussLegendre and NewtonCotes taken from Laurent Fousse's CRQ (Correctly Rounded Quadrature) Numerical Integration code written in MPFR as described in this paper. These algorithm may be selected via the Variant operator as in
⎕FPC←128 {1+1○⍵}∫⍠'g' ○2x ⍝ GaussLegendre 6.28318530717958647692528676655900576839 {1+1○⍵}∫⍠'n' ○2x ⍝ NewtonCotes 6.28318530717958647692528676655900571385 ○2x ⍝ Exact answer 6.28318530717958647692528676655900576839
The Order of the Numerical Integration (the number of rectangles used to approximate the result) is, by default, 128 for the GaussLegendre algorithm, and 32 for NewtonCotes. That number may be changed via the Variant operator as in
N←17x ⎕PP←60 ⋄ ⎕FPC←512 {÷1+⍵*2}∫⍠'g' N 1.51204050407917392632913838918797965662193428158710826434397 {÷1+⍵*2}∫⍠'n' N 1.51204050407917392632913834207540700696171248671126465236073 {÷1+⍵*2}∫⍠('n' 64) N 1.51204050407917392632913838918797965662193428158722152948792 {÷1+⍵*2}∫⍠('n' 96) N 1.51204050407917392632913838918797965662193428158710826434397 ¯3○N ⍝ Exact answer
There are additional choices for algorithms such as
Algorithm  Syntax  Source 

GaussLegendre  ∫⍠'G'  David H. Bailey's ARPREC ARPREC C++ code 
Error Function  ∫⍠'E'  David H. Bailey's ARPREC ARPREC C++ code 
TanhSinh  ∫⍠'T'  David H. Bailey's ARPREC ARPREC C++ code 
TanhSinh  ∫⍠'t' and ∫ as it is the default algorithm 
David H. Bailey's ARPREC ARPREC C++ code translated to C 
GaussLegendre  ∫⍠'g'  Laurent Fousse's CRQ code 
NewtonCotes  ∫⍠'n'  Laurent Fousse's CRQ code 
TanhSinh  ∫⍠'d'  Graeme Dennes Excel Spreadsheet code 
When the Variant operator is applied to the default algorithm, the first number is used to replace the exponent N=5 in the loop limit of 12×2*N.
Not all of these additional algorithm's may be supported in future versions.
Examples
The formula for the Unit Parabola is {⍵*2} whose Integral is {1r3×⍵*3}, and so the Integral of the Unit Parabola from a to a is /{1r3×⍵*3}a,a ←→ 2r3×a*3:
⎕FPC←128 a←2 ⋄ (a){⍵*2}∫a 5.33333333333333333333333333333333333334
This may be used to calculate the Area Enclosed Between a Parabola and a Chord, the formula for which is twothirds of the area of the surrounding rectangle of the chord, the line tangent to the Vertex, and the two sides x=a and x=a.
Because the Unit Parabola curves upward and away from the Xaxis, the above Integral returns the area under (i.e. below and outside) the Parabola. In this case, because we're Integrating from a to a, the chord is the line y={⍵*2} a ←→ y=a*2. The top of the rectangle (chord) is of length aa ←→ 2×a and of height y ←→ a*2, whose area is (2×a)×a*2 ←→ 2×a*3. Subtracting the result of the Integral from the rectangle's area and multiplying this by twothirds yields the area enclosed between a chord and a Parabola.
2r3×(2×a*3)(a){⍵*2}∫a
7.1111111111111111111111111111111111111
The Integral of the 1+Sine function from 0 to ○2 is ○2:
{1+1○⍵}∫○2x ⋄ ○2x 6.28318530717958647692528676655900576834 6.28318530717958647692528676655900576839
and the Integral of the Sine function from 0 to ○2 is (essentially) 0
{1○⍵}∫○2x 6.69453502284549420282962009846618284241E¯42
A Normal Distribution is defined as nd←{(*¯0.5×⍵*2)÷√○2x}. Because this instance of the algorithm is scaled properly, integrating it over the entire width from ¯∞ to ∞ yields an answer of 1 (the area under the curve). However, to achieve this result, we need to increase the default loop limit exponent:
¯∞ nd∫ ∞ 1.000000000000002522695227400560226259006 ¯∞ nd∫⍠10 ∞ 1
Integrating this same function for one, two, and three standard deviations on either side yields the 3sigma rule of :
⍪¯1 ¯2 ¯3 nd∫¨ 1 2 3 0.682689492137085897170465091264075844958 ⍝ 68% 0.954499736103641585599434725666933125059 ⍝ 95% 0.997300203936739810946696370464810045244 ⍝ 99.7%
which describes about how many of the values in a normal distribution lie within one, two, and three standard deviations from the mean.
Numerical Differentiation
Note that Derivative, the inverse of this operator, has also been implemented. These two features may be combined in various ways, for example, to calculate Arc Length.
The formula for Arc Length is simply arclen←{4○⍺⍺∂⍵} which can be used on any given function. For example, the upper half of a unit circle is represented by the function {0○⍵}. The interval between [√0.5, √0.5] is then a quarter of a circle, whose length is
⎕FPC←512 ⋄ ⎕PP←60 (√0.5x){0○⍵}arclen∫√0.5x 1.57079632679489661923132169163975144209858469968755291048747 ○0.5x ⍝ Exact length 1.57079632679489661923132169163975144209858469968755291048747