2018年4月1日 星期日

Numerical Integration


Numerical Integration



Ching-Tang Tseng
Hamilton, New Zealand
2 April 2018
公告:

自2018年元月起,本人拒用臉書,並刪除過去十年內所有貼文。


目前改用嘰嘰喳喳(Twitter)社群網頁,個人帳號使用真名CT Tseng。

40 years ago, our teacher of numerical analysis told us that answer’s precision of numerical integrations were better than numerical derivatives after computer calculation. Those years, main frame computer dominated our digital world. To prove such a calculating precision was not so easy to us. Our computer using times were under controlled rigidly. Each student could use computer once per week only.

Now, I am not only running many computers all day long, but also able to design ABC FORTH functions on many FORTH systems. It is my hobby to see the calculation precision features in my personal ABC FORTH system. In computer calculations, I preferred having good precision rather than to get nice speed.


(1). Introduction


There was a Simpson’s method numerical definite integration example on the Rosetta code website, thanks to it. Even though many textbooks of numerical analysis can supply the Simpson’s method material, for the moment, everybody is still can very convenient to refer to Rosetta code website. We don’t have to so stingy to say thanks to this website when we are citing its contribution. I copy the question and its suggestion directly from that website page as following:



********** Quoted material is beginning here. **********
Question:
Write functions to calculate the definite integral of a function f(x) using all five of the following methods:

rectangular
left
right
midpoint
trapezium
Simpson's

Your functions should take in the upper and lower bounds (a) and (b), and the number of approximations to make in that range (n).

Assume that your example already has a function that gives values for f(x).

Simpson's method is defined by the following pseudo-code:

h := (b - a) / n
sum1 := f(a + h/2)
sum2 := 0

loop on i from 1 to (n - 1)
sum1 := sum1 + f(a + h * i + h/2)
sum2 := sum2 + f(a + h * i)

answer := (h / 6) * (f(a) + f(b) + 4*sum1 + 2*sum2)
Demonstrate your function by showing the results for:

f(x) = x3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25.
f(x) = 1/x, where x is [1,100], with 1,000 approximations. The exact result is the natural log of 100, or about 4.605170
f(x) = x, where x is [0,5000], with 5,000,000 approximations. The exact result is 12,500,000.
f(x) = x, where x is [0,6000], with 6,000,000 approximations. The exact result is 18,000,000.

********** Quoted material is ending here. **********



Four suggested mathematical functions listed above, only this one: f(x) = 1/x is my favorite.
And I would like to add another mathematical function: f(x) = 1/(1+x**2) to be discussed in this article.

(2). Code


My ABC FORTH code is as following:



\ Integrations in lina64 ABC FORTH

\ Integration.f
\ Author: Ching-Tang Tseng
\ 20161107, Hamilton, NZ

9 reals a b h sum1 sum2 x f(x) f(a) f(b)
2 integers n i
1000 value approximations

: 1/x ( -- )
basic
10 let { f(x) = 1. E 0 / x }
20 end ;

: 1/(1+x**2) ( -- )
basic
10 let { f(x) = 1. E 0 / ( 1.0 e 0 + x * x ) }
20 end ;
defer function
' 1/x is function

: range1/x ( -- )
basic
10 let n = approximations
20 let { a = f1.0e0 }
30 let { b = 100.0 e 0 }
40 end ;

: range1/(1+x**2) ( -- )
basic
10 let n = approximations
20 let { a = f0.0e0 }
30 let { b = f2.0e0 }
40 end ;
defer setting
' range1/x is setting
\ All parts above are changeable setting.

\ All parts below are fixed.
: SimpsonInit ( -- )
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }

30 let { x = a + h / 2. E 0 }
40 run function
50 let { sum1 = f(x) }

60 let { sum2 = 0. E 0 }

70 let { x = a }
80 run function
90 let { f(a) = f(x) }

100 let { x = b }
110 run function
120 let { f(b) = f(x) }

130 end ;

: Simpson ( -- )
basic
10 run SimpsonInit

20 for i = 1 to n - 1

30 let { x = a + h * I>R ( i ) + h / 2. E 0 }
40 run function
50 let { sum1 = sum1 + f(x) }

60 let { x = a + h * I>R ( i ) }
70 run function
80 let { sum2 = sum2 + f(x) }

90 next i

100 print " Simpson integration = " ;
{ ( h / 6. E 0 ) * ( f(a) + f(b) + 4. E 0 * sum1 + 2. E 0 * sum2 ) }

110 end ;

: L-Rectangle ( -- )
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a }
40 let { sum1 = 0. E 0 }

50 for i = 1 to n
60 run function
70 let { sum1 = sum1 + h * f(x) }
80 let { x = x + h }
90 next i

100 print " L-Rectangle integration = " ; { sum1 }

110 end ;

: R-Rectangle ( -- )
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a }
40 let { sum1 = 0. E 0 }

50 for i = 1 to n
60 let { x = x + h }
70 run function
80 let { sum1 = sum1 + h * f(x) }
90 next i

100 print " R-Rectangle integration = " ; { sum1 }

110 end ;

: M-Rectangle ( -- )
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }
30 let { x = a }
40 let { sum1 = 0. E 0 }

50 for i = 1 to n
60 let { x = x + h / 2. E 0 }
70 run function
80 let { sum1 = sum1 + h * f(x) }
90 let { x = x + h / 2. E 0 }
100 next i

110 print " M-Rectangle integration = " ; { sum1 }

120 end ;

: Trapezoidal ( -- )
basic
10 run Setting

20 let { h = ( b - a ) / I>R ( n ) }
30 let { sum1 = 0. E 0 }

40 let { x = b }
50 run function
60 let { f(b) = f(x) }
70 let { x = a }
80 run function
90 let { f(a) = f(x) }
100 let { sum1 = h * ( f(b) + f(a) ) / 2. E 0 }

110 for i = 1 to n - 1
120 let { x = x + h }
130 run function
140 let { sum1 = sum1 + h * f(x) }
150 next i

160 print " Trapezoidal integration = " ; { sum1 }

170 end ;

\ Output description depends on input setting.

: integrations ( -- )
L-Rectangle
R-Rectangle
M-Rectangle
Trapezoidal
Simpson
;

: function1
cr cr
." Integration of 1/x from 1 to 100 = ln(100) "
cr
integrations
cr cr 18 spaces
." ln(100) = " {{ ln ( 100.0 e 0 ) }}
fs. cr ;

: function2 ( -- )
cr cr
." Integration of 1/(1+x^2) from 0 to 2 = atan(2) "
cr
integrations
cr cr 18 spaces
." atan(2) = " {{ atan ( f2.0e0 ) }}
fs. cr ;

: (main)
['] 1/x is function ['] range1/x is setting function1
['] 1/(1+x**2) is function ['] range1/(1+x**2) is setting function2
;

: main
cr ." =====Integration===== "
cr ." approximations = 1000 "
1000 to approximations (main)
cr ." ========================================= "
cr ." approximations = 10000 "
10000 to approximations (main)
cr ." =======the end======= "
;

main


(3). Execution Outcomes



=====Integration=====
approximations = 1000

Integration of 1/x from 1 to 100 = ln(100)

L-Rectangle integration = 4.654991057514675655 E 0
R-Rectangle integration = 4.556981057514675656 E 0
M-Rectangle integration = 4.604762548678374701 E 0
Trapezoidal integration = 4.605986057514675655 E 0
Simpson integration = 4.605170384957141681 E 0

ln(100) = 4.605170185988091362 E 0


Integration of 1/(1+x^2) from 0 to 2 = atan(2)

L-Rectangle integration = 1.107948664460762084 E 0
R-Rectangle integration = 1.106348664460762086 E 0
M-Rectangle integration = 1.107148744460752491 E 0
Trapezoidal integration = 1.107148664460762086 E 0
Simpson integration = 1.107148717794089129 E 0

atan(2) = 1.107148717794090466 E 0

=========================================
approximations = 10000

Integration of 1/x from 1 to 100 = ln(100)

L-Rectangle integration = 4.610078852591290491 E 0
R-Rectangle integration = 4.600277852591290491 E 0
M-Rectangle integration = 4.60516610271650127 E 0
Trapezoidal integration = 4.605178352591290491 E 0
Simpson integration = 4.605170186008097743 E 0

ln(100) = 4.605170185988091362 E 0


Integration of 1/(1+x^2) from 0 to 2 = atan(2)

L-Rectangle integration = 1.107228717260755149 E 0
R-Rectangle integration = 1.107068717260755151 E 0
M-Rectangle integration = 1.107148718060755166 E 0
Trapezoidal integration = 1.10714871726075515 E 0
Simpson integration = 1.107148717794089581 E 0

atan(2) = 1.107148717794090466 E 0

=======the end======= OK



(4). Discussion




We are familiar with the following integral formulas:

1/x dx = ln(x)
1/(1+x**2) dx = arctan(x)

In general, functions ln(x) and arctan(x) are always to be designed as primitive functions. So we are able to check the precision between the numerical integration calculation outcomes and the system function values itself.

Integral techniques are still having many other methods. Here we have done Simpson’s method and some other very fundamental methods. All output results are comparable.

I am using Lina64 under Ubuntu 16.04 OS. Emphasized once again, pure software floating point system was designed by me. Only 15 digits precision are reliable, but I am still printing out true values with total digits frankly.

What our teacher had ever told us about the precision of numerical integration outcomes from the computer calculation? After 40 years later, we are able to look at the real calculation results via FORTH quickly and easily.

Seeing that precision of numerical integration could reach to a high precision level, can we have it to be a method to supply us both FLN and FATAN functions? According to the fundamental theorem of integral calculus, the answer is affirmative yes. Problem is: How many approximations must be taken? This article shows us that at least few thousand times must to take. But in general, to get a primitive function value, carry out by Taylor series method, usually needs only less than one hundred times.

Lina64 is powerful enough. There were many features been demonstrated in the past articles already. This time, I would like to show you Lina64’s graphic ability in video. I did not design graphic package in Lina64. Other people had contributed a lot to Ubuntu already. Lina64 is able to manipulate all kinds of interpreter style language, for example, Python, Perl...etc., Graphic package included as well. This convenience makes me have the opportunity to show out quickly many mathematical function figures on the screen. Those graphic codes are out of FORTH code discussing range. Their coding format are the same as in the Win32Forth system, to be treated as an Object-Oriented programming.

Graphic display is helpful to explain mathematical problem. I would like to post 4 figures of function: f(x)=1/(1+x**2), f(x)=arctan(x), f(x)=1/x, f(x)=ln(x). They are related to the topic talked above in this article. The quality of a mathematical diagram is highly depended on the spanned range you selected. To adjust out a suitable setting have to spend plenty of time. This video shows that Lina64 can help a lot.



Fig.1





Fig.2






Fig.3



Fig.4





(5). Demonstrated Video





(6). My personal testing remarks

primitive functions FLN, FATAN approved.






沒有留言: