Numerical Integration
Ching-Tang
Tseng
Hamilton,
New Zealand
2
April 2018
公告：
自2018年元月起，本人拒用臉書，並刪除過去十年內所有貼文。
目前改用嘰嘰喳喳(Twitter)社群網頁，個人帳號使用真名CT Tseng。
自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
沒有留言:
張貼留言