118 lines
3.6 KiB
Plaintext
118 lines
3.6 KiB
Plaintext
while 1
|
|
read x$
|
|
if x$ ="end" then print "**Over**": end
|
|
|
|
read a, b, N, knownValue
|
|
|
|
print " Function y ="; x$; " from "; a; " to "; b; " in "; N; " steps"
|
|
print " Known exact value ="; knownValue
|
|
|
|
areaLR = IntegralByLeftRectangle( x$, a, b, N)
|
|
areaRR = IntegralByRightRectangle( x$, a, b, N)
|
|
areaMR = IntegralByMiddleRectangle( x$, a, b, N)
|
|
areaTr = IntegralByTrapezium( x$, a, b, N)
|
|
areaSi = IntegralBySimpsonRule( x$, a, b, N)
|
|
|
|
print "Left rectangle method "; using( "##########.##########", areaLR); " diff "; knownValue-areaLR; tab(70); (knownValue-areaLR)/knownValue*100;" %"
|
|
print "Right rectangle method "; using( "##########.##########", areaRR); " diff "; knownValue-areaRR; tab(70); (knownValue-areaRR)/knownValue*100;" %"
|
|
print "Middle rectangle method "; using( "##########.##########", areaMR); " diff "; knownValue-areaMR; tab(70); (knownValue-areaMR)/knownValue*100;" %"
|
|
print "Trapezium method "; using( "##########.##########", areaTr); " diff "; knownValue-areaTr; tab(70); (knownValue-areaTr)/knownValue*100;" %"
|
|
print "Simpson's Rule "; using( "##########.##########", areaSi); " diff "; knownValue-areaSi; tab(70); (knownValue-areaSi)/knownValue*100;" %"
|
|
|
|
print
|
|
|
|
wend
|
|
|
|
end
|
|
|
|
'------------------------------------------------------
|
|
'we have N sizes, that gives us N+1 points
|
|
'point 0 is a
|
|
'point N is b
|
|
'point i is xi =a +i *h
|
|
'Often, precision is (sharper?) then single step area
|
|
'So there should be EXACT number of steps, hence loop by integer i.
|
|
|
|
function IntegralByLeftRectangle( x$, a, b, N)
|
|
h = ( b -a) /N
|
|
s = 0
|
|
for i = 0 to N -1
|
|
x = a +i *h
|
|
s = s + h *eval( x$)
|
|
next
|
|
IntegralByLeftRectangle = s
|
|
end function
|
|
|
|
function IntegralByRightRectangle( x$, a, b, N)
|
|
h =( b -a) /N
|
|
s = 0
|
|
for i =1 to N
|
|
x = a +i *h
|
|
s = s + h *eval( x$)
|
|
next
|
|
IntegralByRightRectangle = s
|
|
end function
|
|
|
|
function IntegralByMiddleRectangle( x$, a, b, N)
|
|
h =( b -a) /N
|
|
s = 0
|
|
for i =0 to N -1
|
|
x = a +i *h +h /2
|
|
s = s + h *eval( x$)
|
|
next
|
|
IntegralByMiddleRectangle = s
|
|
end function
|
|
|
|
function IntegralByTrapezium( x$, a, b, N)
|
|
'Formula is h*((f(a)+f(b))/2 + sum_{i=1}^{N-1} (f(x_i)))
|
|
h =( b -a) /N
|
|
x = a
|
|
fa =eval( x$)
|
|
x =b
|
|
fb =eval( x$)
|
|
s = h *( fa +fb) /2
|
|
for i =1 to N -1
|
|
x = a +i *h
|
|
s = s + h *eval( x$)
|
|
next
|
|
IntegralByTrapezium = s
|
|
end function
|
|
|
|
function IntegralBySimpsonRule( x$, a, b, N)
|
|
'Simpson
|
|
'N should be even.
|
|
if N mod 2 then N =N +1
|
|
'It really doesn't look right to double number of points from N to 2N -
|
|
' - this method is most accurate of all presented!
|
|
'So we use NN as N/2, and N will be 2NN
|
|
'Formula is h/6*( f(a)+f(b) + 4*(f(x_1)+f(x_3)+...+f(x_{2NN-1})+ 2*(f(x_2)+f(x_4)+...+f(x_{2NN-2})) )
|
|
'Somehow I messed up h/6, h/3 and what is h, regarding "n=number of double intervals of size 2h"
|
|
NN =N /2
|
|
|
|
h =( b -a) /N
|
|
x =a
|
|
fa =eval (x$)
|
|
x =b
|
|
fb =eval( x$)
|
|
s = h /3 *( fa +fb)
|
|
for i =1 to 2 *NN -1 step 2
|
|
x = a +i *h
|
|
s = s + h /3 *4 *eval( x$) 'odd points
|
|
next
|
|
for i =2 to 2 *NN -2 step 2
|
|
x = a +i *h
|
|
s = s + h /3 *2 *eval( x$) 'even points
|
|
next
|
|
|
|
IntegralBySimpsonRule = s
|
|
end function
|
|
|
|
'=======================================================
|
|
data "x^3", 0, 1, 100, 0.25
|
|
data "x^-1", 1, 100, 1000, 4.605170
|
|
data "x", 0, 5000, 1000, 12500000.0 ' should use 5 000 000 steps
|
|
data "x", 0, 6000, 1000, 18000000.0 ' should use 6 000 000 steps
|
|
data "end"
|
|
|
|
end
|