一百個例題 (71 ~ 75)
Ching-Tang Tseng
Hamilton, New Zealand
28 September 2024
\ Dmultiply.f \ In FSL_UTIL the following is defined: : d*1 ( d1 d2 -- dprod ) \ double multiply dup 3 pick xor >r dabs 2swap dabs udm* r> 0= IF exit ENDIF dxor 2swap dxor 1. UD+c >r \ ud1 ud2 carry1 -- ud1+ud2 carry2 2swap r> UD+c drop ; : d*2 ( d1 d2 -- dprod ) >r >r over r@ um* rot r> * + rot r> * + ; : d*3 ( d1 d2 -- d1*d2 ) >r swap >r 2dup um* 2swap r> * swap r> * + + ; : d*4 ( multiplicand . multiplier . -- product . ) 3 PICK * >R TUCK * >R UM* R> + R> + ;
\ FILEOP.F \ This file operation program works in Lina64 under CentOS7 \ Author: Ching Tang Tseng \ Date : 20160310, Hamilton, New Zealand WANT ALLOCATE \ -rw0rw0rw0 = 110110110b = 438d : all are R/W enable 438 CONSTANT R/W VARIABLE FileID VARIABLE Fptr VARIABLE Frem VARIABLE Flen VARIABLE Fsize 10240 Fsize ! \ (1)floating Fadr \ : Fadr PAD 4096 + ; \ (2)allocate Fadr ??? \ Fsize ALLOCATE DROP CONSTANT Fadr \ (3)fixed Fadr: 100 KB below EM \ EM HERE - . --> get 33425420 --> 33 MB free spaces : Fadr EM 102400 - ; : SetUpFptrFrem ( -- ) Fadr Fptr ! Flen @ Frem ! ; : (FILE.) ( -- addr len ) Fadr Flen @ ; : FileType ( -- ) CR (FILE.) TYPE ; : FileDump ( -- ) CR (FILE.) DUMP ; \ Beware! only a R/W attributed file can be manipulated \ S" Filename.f" GetFile = "Filename.f" GetFile : GetFile ( addr len -- ) Fadr Fsize @ 0 FILL R/W OPEN-FILE IF CR ABORT" OPEN-FILE error?" THEN FileID ! CR ." File ID is : " FileID @ . Fadr Fsize @ FileID @ READ-FILE IF CR ABORT" READ-FILE error?" THEN DUP Flen ! CR . ." Bytes has been read!" FileID @ CLOSE-FILE IF CR ABORT" CLOSE-FILE error!" THEN SetUpFptrFrem ; \ use PAD area create all blanks : NewFile ( addr len n -- ) Flen ! PAD Flen @ 32 FILL R/W CREATE-FILE IF CR ABORT" CREATE-FILE error!" THEN FileID ! PAD Flen @ FileID @ WRITE-FILE IF CR ABORT" WRITE-FILE error!" THEN CR Flen @ . ." Bytes has been written!" FileID @ CLOSE-FILE IF CR ABORT" CLOSE-FILE error!" THEN SetUpFptrFrem ; \ Beware! Flen must be set, before you WriteFile : SaveFile ( addr len -- ) R/W CREATE-FILE IF CR ABORT" CREATE-FILE error!" THEN FileID ! CR ." FileID is: " FileID @ . Fadr Flen @ FileID @ WRITE-FILE IF CR ABORT" WRITE-FILE error!" THEN CR Flen @ . ." Bytes has been written!" FileID @ CLOSE-FILE IF CR ABORT" CLOSE-FILE error!" THEN ; \ for simple testing usage: \ S" This is a simple test." SendText>F \ : SendText>F ( adr n -- ) \ DUP Flen ! Fadr SWAP MOVE ; \ Frem, Fptr are to be used for other testing.
\ PRINTstudy.f
: main1 ( -- )
cr
73 65
do
13 1
do
j emit i .
loop
cr
loop ;
2 integers i j
: main2 ( -- )
basic
10 run cr
20 for j = 65 to 72
30 for i = 1 to 12
40 run j emit i .
50 next i
60 run cr
70 next j
80 end ;
\s
main1
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12
B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12
D1 D2 D3 D4 D5 D6 D7 D8 D9 D10 D11 D12
E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 E11 E12
F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12
G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12
H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 H11 H12
ok
main2
A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12
B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12
D1 D2 D3 D4 D5 D6 D7 D8 D9 D10 D11 D12
E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 E11 E12
F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12
G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 G11 G12
H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 H11 H12
ok
\ (74)SimpsonIntegration.f
fvariable step
defer method ( fn F: x -- fn[x] )
: left execute ;
: right step f@ f+ execute ;
: mid step f@ 2e f/ f+ execute ;
: trap
dup fdup left
fswap right f+ 2e f/ ;
: simpson
dup fdup left
dup fover mid 4e f* f+
fswap right f+ 6e f/ ;
: set-step ( n F: a b -- n F: a )
fover f- dup 0 d>f f/ step f! ;
: integrate ( xt n F: a b -- F: sigma )
set-step
0e
0 do
dup fover method f+
fswap step f@ f+ fswap
loop
drop fnip
step f@ f* ;
\ testing similar to the D example
: test
' is method ' 4 -1e 2e integrate f. ;
: fn1 fsincos f+ ;
: fn2 fdup f* 4e f* 1e f+ 2e fswap f/ ;
7 set-precision
test left fn2 \ 2.456897
test right fn2 \ 2.245132
test mid fn2 \ 2.496091
test trap fn2 \ 2.351014
test simpson fn2 \ 2.447732
\S
Numerical integration
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.
\ (74-1)RectangleIntegration.f
\ Author: Ching-Tang Tseng
\ 20161107, Hamilton, Hamilton, NZ
9 reals a b h sum1 sum2 x f(x) f(a) f(b)
2 integers n i
: function ( -- )
basic
10 let { f(x) = 1 / x }
20 end ;
: Setting ( -- )
basic
10 let n = 1000
20 let { a = 1 }
30 let { b = 100 }
40 end ;
\ All parts above are adjustable 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 }
40 run function
50 let { sum1 = f(x) }
60 let { sum2 = 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 }
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 ) * ( f(a) + f(b) + 4 * sum1 + 2 * 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 }
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 }
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 }
50 for i = 1 to n
60 let { x = x + h / 2 }
70 run function
80 let { sum1 = sum1 + h * f(x) }
90 let { x = x + h / 2 }
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 }
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 }
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 ;
18 sigdigits !
: main ( -- )
cr cr ." Integration of 1/x from 1 to 100 = Ln(100) " cr
L-Rectangle
R-Rectangle
M-Rectangle
Trapezoidal
Simpson
cr cr ." Ln(100) = " 100.0e0 fln f. cr ;
main
\S
Integration of 1/x from 1 to 100 = Ln(100)
L-Rectangle integration = 4.65499105751467615
R-Rectangle integration = 4.55698105751467615
M-Rectangle integration = 4.60476254867837518
Trapezoidal integration = 4.60598605751467615
Simpson integration = 4.60517038495714217
Ln(100) = 4.60517018598809137
ok
\ (75)高斯消去法解多元聯立一次方程式2009-08-22
INTEGER #ROW INTEGER ITEMP \ #ROW為方陣量,本例為3X3矩陣
INTEGER M INTEGER M1 INTEGER M2
REAL FTEMP
10 ARRAY Y 10 ARRAY X \ 矩陣上限宣告,暫設為10X10
10 10 MATRIX AA
: 輸入數據
[[ #ROW = 3 ]]
{{ AA ( 1 1 ) = 3 }} {{ AA ( 1 2 ) = 2 }} {{ AA ( 1 3 ) = 1 }}
{{ AA ( 2 1 ) = 2 }} {{ AA ( 2 2 ) = 4 }} {{ AA ( 2 3 ) = 5 }}
{{ AA ( 3 1 ) = 4 }} {{ AA ( 3 2 ) = 1 }} {{ AA ( 3 3 ) = 1 }}
{{ Y ( 1 ) = 17 }} {{ Y ( 2 ) = 41 }} {{ Y ( 3 ) = 16 }} ;
: 高斯處理
#ROW 1
DO
[[ M = I + 1 ]]
#ROW 1+ M
DO
#ROW 1+ M
DO
{{ AA ( J I ) = AA ( J I ) - ( AA ( K I ) * AA ( J K ) ) / AA ( K K ) }}
LOOP
{{ Y ( I ) = Y ( I ) - AA ( I J ) * Y ( J ) / AA ( J J ) }}
LOOP
M #ROW =
IF {{ FTEMP = AA ( M M ) }}
FTEMP F0=
IF ." Equation disparity! " QUIT THEN
THEN
LOOP ;
: 得到解答
[[ M = #ROW + 1 ]]
#ROW 1+ 1
DO
[[ M1 = M - I ]]
{{ X ( M1 ) = Y ( M1 ) / AA ( M1 M1 ) }}
[[ ITEMP = ( M - I ) - #ROW ]]
ITEMP DUP
0>
IF ." Trival! " QUIT THEN
0<
IF #ROW 1+ 1
DO
[[ M2 = M - I ]]
{{ X ( M1 ) = X ( M1 ) - AA ( M1 M2 ) * X ( M2 ) / AA ( M1 M1 ) }}
[[ ITEMP = ( M - I - 1 ) - ( M - J ) ]]
ITEMP
0=
IF LEAVE THEN
LOOP
THEN
LOOP ;
: 輸出答案
#ROW 1+ 1
DO
{{ FTEMP = X ( I ) }}
CR ." X( " I . ." )= " FTEMP F.
LOOP ;
: GAUSS
輸入數據
高斯處理
得到解答
輸出答案
;
沒有留言:
張貼留言