一百個例題 (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 輸入數據 高斯處理 得到解答 輸出答案 ;
沒有留言:
張貼留言