一百個例題 (37 ~ 40)
Ching-Tang Tseng
Hamilton, New Zealand
5 September 2024
FORTH:慶的世界。 本網頁為中華民國(Republic Of China)FORTH語言協會 永久性網頁
2024年9月4日 星期三
第(37)個範例，是一個完全在複數體系內運作的展示程式。
我在網上公佈我的系統具有獨立複數體系的程式運作功能前，全世界還找不到這方面的公開資料，也沒有任何程式語言可以寫這種程式來進行測試。以前，大家在使用電腦求解複數體系的問題時，都使用只具有浮點實數計算功能的系統設計程式。現在，可能很多程式語言都已具備有複數體系的程式功能了。我未作研究，不作評論。
從唸初中開始，我們的數學課本就告訴我們：
一元二次方程式的標準通式是： aX^2+bX+c=0
其兩個根的標準公式就是： 『2a分之負b加減根號b平方減4ac』
因為這個公式很簡單，大家也容易了解，所以我經常使用它來展示系統性能，這個展示複數體系的範例程式也不例外。
問題看似非常簡單，不實作的人，就很難體會其中的奧妙。一個初中生就能懂得的問題，其道理卻蘊含了數學體系性質的延伸性，我曾經問過自己：在浮點數體系中，沒有疑問的公式，能不能照樣運用在複數體系中？這樣的問題，就是我所謂之數學體系性質的延伸性。沒有實作出結果前，不要亂下結論，我做出來了，就敢直說：可以。
在學時，沒有老師敢告訴大家這種性質，因為，想測試複數體系的數學計算，非常冗長又麻煩，在有限的教學時間內，老師只會教升學考試要考的東西，理論性的材料，大部分都只點到為止的講一講。前面第(32)個範例，講過算出複數的平方根的程序很麻煩，光是只取主值來代表就很難算了，老師怎麼可能用冗長的計算來煩學生也煩到了自己？
這並不是只有華人的環境才如此，因為我有個 kForth 系統的作者朋友，他是印度裔的加拿大高能物理博士，跟我書信來往時，發現了我提到的這種問題，他就純用 Forth 寫了一個程式，捐給 FSL(Forth Scientifical Library) 當庫存公益程式， FSL 組織也接受了。這就表示，全世界以前確實是沒有人注意到這種問題，要等到電腦系統很容易使用後，才能以數值解的方式把問題表達清楚。
這個範例程式很簡單，所以內容請自己看。但是，請注意範例執行中的過程與結果。
首先，請注意複數計算環境，是被兩個中括號 『 [ ----- ] 』包裹起來的計算環境。
其次，過程中告訴大家，複數在透明的 Forth 系統中，真實的結構只是兩個浮點數，而且不帶有那個表示虛數的 i 。電腦系統只在需要印出結果給人類看時，才尾綴一個 i 來表示它是個複數，一般程式語言很難設計這種格式的輸入方法，所以都把 i 放在第二個浮點數前面來處理數字，我設計的系統就將 i 放在課堂上與書本中同樣的標準位置來處理。所以，您在對答式輸入的過程中，每次只輸入兩個浮點數，而我又在 ABC Forth 系統中增添了設計，在使用者輸入整數或浮點數進入複數系統時，會自動全部轉換成純粹的複數，這方面的設計技術相當複雜，我不擬在此討論。
看完了程式，也請執行出結果，您就會注意到，我寫了回算的驗證程式，結果竟然不是 0 。不要感到奇怪，我的設計也沒有錯，答案就是這樣。只需注意核算結果所得數字，是不是小到 10 的負十幾次方去了？是，那就表示是 0 。這是電腦處理浮點數字的固有本質，更何況，無數的方程式，其根都是無法以有限位數表達出來的無理數，可能連這個範例中的解答也是，但我沒用我的大數系統去驗證，我可以驗證得出來。
這麼簡單的一個連初中生都學過的問題，怎麼牽扯出了這許多的討論？答案是我重視實作，實踐才是檢驗真理的唯一標準，確實如此。
:
我在網上公佈我的系統具有獨立複數體系的程式運作功能前，全世界還找不到這方面的公開資料，也沒有任何程式語言可以寫這種程式來進行測試。以前，大家在使用電腦求解複數體系的問題時，都使用只具有浮點實數計算功能的系統設計程式。現在，可能很多程式語言都已具備有複數體系的程式功能了。我未作研究，不作評論。
從唸初中開始，我們的數學課本就告訴我們：
一元二次方程式的標準通式是： aX^2+bX+c=0
其兩個根的標準公式就是： 『2a分之負b加減根號b平方減4ac』
因為這個公式很簡單，大家也容易了解，所以我經常使用它來展示系統性能，這個展示複數體系的範例程式也不例外。
問題看似非常簡單，不實作的人，就很難體會其中的奧妙。一個初中生就能懂得的問題，其道理卻蘊含了數學體系性質的延伸性，我曾經問過自己：在浮點數體系中，沒有疑問的公式，能不能照樣運用在複數體系中？這樣的問題，就是我所謂之數學體系性質的延伸性。沒有實作出結果前，不要亂下結論，我做出來了，就敢直說：可以。
在學時，沒有老師敢告訴大家這種性質，因為，想測試複數體系的數學計算，非常冗長又麻煩，在有限的教學時間內，老師只會教升學考試要考的東西，理論性的材料，大部分都只點到為止的講一講。前面第(32)個範例，講過算出複數的平方根的程序很麻煩，光是只取主值來代表就很難算了，老師怎麼可能用冗長的計算來煩學生也煩到了自己？
這並不是只有華人的環境才如此，因為我有個 kForth 系統的作者朋友，他是印度裔的加拿大高能物理博士，跟我書信來往時，發現了我提到的這種問題，他就純用 Forth 寫了一個程式，捐給 FSL(Forth Scientifical Library) 當庫存公益程式， FSL 組織也接受了。這就表示，全世界以前確實是沒有人注意到這種問題，要等到電腦系統很容易使用後，才能以數值解的方式把問題表達清楚。
這個範例程式很簡單，所以內容請自己看。但是，請注意範例執行中的過程與結果。
首先，請注意複數計算環境，是被兩個中括號 『 [ ----- ] 』包裹起來的計算環境。
其次，過程中告訴大家，複數在透明的 Forth 系統中，真實的結構只是兩個浮點數，而且不帶有那個表示虛數的 i 。電腦系統只在需要印出結果給人類看時，才尾綴一個 i 來表示它是個複數，一般程式語言很難設計這種格式的輸入方法，所以都把 i 放在第二個浮點數前面來處理數字，我設計的系統就將 i 放在課堂上與書本中同樣的標準位置來處理。所以，您在對答式輸入的過程中，每次只輸入兩個浮點數，而我又在 ABC Forth 系統中增添了設計，在使用者輸入整數或浮點數進入複數系統時，會自動全部轉換成純粹的複數，這方面的設計技術相當複雜，我不擬在此討論。
看完了程式，也請執行出結果，您就會注意到，我寫了回算的驗證程式，結果竟然不是 0 。不要感到奇怪，我的設計也沒有錯，答案就是這樣。只需注意核算結果所得數字，是不是小到 10 的負十幾次方去了？是，那就表示是 0 。這是電腦處理浮點數字的固有本質，更何況，無數的方程式，其根都是無法以有限位數表達出來的無理數，可能連這個範例中的解答也是，但我沒用我的大數系統去驗證，我可以驗證得出來。
這麼簡單的一個連初中生都學過的問題，怎麼牽扯出了這許多的討論？答案是我重視實作，實踐才是檢驗真理的唯一標準，確實如此。
:
\ (37)複數係數一元二次方程式的根.f
\ Two roots of complex coefficient quadratic equation
\ 2011-04-12 Ching-Tang Tseng, Hamilton NZ
\ Equation: F(x)=A*x^2+B*x+C=0 roots: x1,x2
\ Renww date: 2014-08-08
8 complexs A B C Disc x1 x2 F(x1) F(x2)
: Reset ( -- )
BASIC
10 let [ A = 0 + 0 i ] :: [ B = 0 + 0 i ] :: [ C = 0 + 0 i ]
:: [ DISC = 0 + 0 i ] :: [ X1 = 0 + 0 i ] :: [ X2 = 0 + 0 i ]
20 end
;
: QUAD ( -- )
BASIC
10 let [ Disc = SQRT ( B * B - ( 4 + 0 i ) * A * C ) ]
:: [ x1 = ( ( NEGATE ( B ) + Disc ) ) / ( ( 2 + 0 i ) * A ) ]
:: [ x2 = ( ( NEGATE ( B ) - Disc ) ) / ( ( 2 + 0 i ) * A ) ]
:: [ F(x1) = A * x1 * x1 + B * x1 + C ]
:: [ F(x2) = A * x2 * x2 + B * x2 + C ]
20 end
;
: EnterCoefficients ( -- )
BASIC
10 print " Enter coefficient A, ( 3 - 5 i ) must type: 3 -5 "
20 InputZ A
30 run CR
40 print " Enter coefficient B, ( -2.5 + 0 i ) must type: -2.5 0 "
50 InputZ B
60 run CR
70 print " Enter coefficient C, ( 0 + 6.3 i ) must type: 0 6.3 "
80 InputZ C
90 run CR
100 end
;
: Main ( -- )
BASIC
10 run CR
20 Print " Quadratic Eq.: A*x^2+B*x+C=0 "
30 run CR
=> Reset
=> EnterCoefficients
=> QUAD
=> CR
40 print [ " x1 = " ; x1 ]
50 print [ " F(x1)= " ; F(x1) ]
60 run CR
70 print [ " x2 = " ; x2 ]
80 print [ " F(x2)= " ; F(x2) ]
90 run CR CR
100 end
;
\s
main
Quadratic Eq.: A*x^2+B*x+C=0
Enter coefficient A, ( 3 - 5 i ) must type: 3 -5 ? 3 -5
Enter coefficient B, ( -2.5 + 0 i ) must type: -2.5 0 ? -2.5 0
Enter coefficient C, ( 0 + 6.3 i ) must type: 0 6.3 ? 0 6.3
x1 = -.876151 + .445031 i
F(x1)= -4.44089E-16 + 3.55271E-15 i
x2 = 1.09674 - .077384 i
F(x2)= -1.33227E-15 + 4.44089E-15 i
第(38)個範例是求解反矩陣(Inverse Matrix)的程式。
在數學計算系統中，矩陣是一個獨立的體系，它的運算也自成一格，但運算意義仍依循由整數、浮點數、複數、發展到矩陣的一貫道統，所有的運算指令都得單獨設計，才能形成完整的體系。矩陣體系內有自己的四則運算，卻不足以完成整個體系演算時的需求，矩陣在用來解決數學問題時，需要一些功能函數協助運作，問題才容易獲得解答，功能函數非常多，求解矩陣的反元素就是其中之一，它的性質就像求得 x 的反元素 1/x 一樣，在計算中經常要用它來簡化演算。
這個範例程式就展示一種求得反矩陣的方法，方法有很多種，這裡採用的是高斯消去法(Gauss-Jordan Elimination)，程式中詳列了資料來源，原始資料中也只提供程式算法不講詳細的理論，參考文獻中可以找到理論根據，原始資料說不在話下，我就也不在話下，有興趣的話，可以單憑高斯消去法術語從網上查到資料。 4 X 4 以下的方矩陣還有固定的數學公式可以直接算出反矩陣各個元素的值來， 4 X 4 以上的方矩陣就無法以簡單的數學公式直接計算出各個元素之值了。這個範例程式則展示了 5 X 5 與 8 X 8 的兩個方矩陣之演算結果，同時程式最後也要核算出原矩陣乘上它的反矩陣必須是對應的單位矩陣這樣的結果 ： [A] * 1/[A] = [1]。程式中為了包括這樣的核驗設計，用到矩陣乘法(Matrix*)此一獨立運算指令，指令的規格也適用於任何維度。
設計矩陣體系的重點在矩陣資料結構的規劃，根據傳統程式語言對矩陣問題的處裡慣例，所有的設計也都只作到能取、存矩陣中的元素來進行數學計算就算了，這一點，我設計的系統已經可以辦到。但在我心目中的理想設計，不應僅只是這樣，我認為矩陣的運算、矩陣的函數全都可以設計成通用指令。例如：這個範例程式最終可以設計成已知一個 [A] 矩陣時，只要執行一個名為 1/[A] 的矩陣函數指令，就能得到他的反矩陣。而且，這樣的矩陣函數指令必須適用於任何維度的矩陣。系統中的矩陣體系能設計成這樣才接近於理想。要想辦到此事，矩陣的資料結構必須規劃成維度必須存在於資料結構中，求 1/[A] 時，此函數指令才能自動到資料結構中取得維度數據，然後才處理出結果。原本在 Win32Forth 中發展出來的 ABC Forth 系統並沒有這樣的規劃，2012 年以後，我轉往 Linux 作業系統的環境發展系統，才開始規劃出這樣的矩陣資料結構。這樣做，在想要印出矩陣內容時也會更為方便，只要能取得維度數據，印出矩陣的執行迴路指標就可以被自動設定，一個印出矩陣的指令就可以適用於印出任何維度的矩陣。
:
在數學計算系統中，矩陣是一個獨立的體系，它的運算也自成一格，但運算意義仍依循由整數、浮點數、複數、發展到矩陣的一貫道統，所有的運算指令都得單獨設計，才能形成完整的體系。矩陣體系內有自己的四則運算，卻不足以完成整個體系演算時的需求，矩陣在用來解決數學問題時，需要一些功能函數協助運作，問題才容易獲得解答，功能函數非常多，求解矩陣的反元素就是其中之一，它的性質就像求得 x 的反元素 1/x 一樣，在計算中經常要用它來簡化演算。
這個範例程式就展示一種求得反矩陣的方法，方法有很多種，這裡採用的是高斯消去法(Gauss-Jordan Elimination)，程式中詳列了資料來源，原始資料中也只提供程式算法不講詳細的理論，參考文獻中可以找到理論根據，原始資料說不在話下，我就也不在話下，有興趣的話，可以單憑高斯消去法術語從網上查到資料。 4 X 4 以下的方矩陣還有固定的數學公式可以直接算出反矩陣各個元素的值來， 4 X 4 以上的方矩陣就無法以簡單的數學公式直接計算出各個元素之值了。這個範例程式則展示了 5 X 5 與 8 X 8 的兩個方矩陣之演算結果，同時程式最後也要核算出原矩陣乘上它的反矩陣必須是對應的單位矩陣這樣的結果 ： [A] * 1/[A] = [1]。程式中為了包括這樣的核驗設計，用到矩陣乘法(Matrix*)此一獨立運算指令，指令的規格也適用於任何維度。
設計矩陣體系的重點在矩陣資料結構的規劃，根據傳統程式語言對矩陣問題的處裡慣例，所有的設計也都只作到能取、存矩陣中的元素來進行數學計算就算了，這一點，我設計的系統已經可以辦到。但在我心目中的理想設計，不應僅只是這樣，我認為矩陣的運算、矩陣的函數全都可以設計成通用指令。例如：這個範例程式最終可以設計成已知一個 [A] 矩陣時，只要執行一個名為 1/[A] 的矩陣函數指令，就能得到他的反矩陣。而且，這樣的矩陣函數指令必須適用於任何維度的矩陣。系統中的矩陣體系能設計成這樣才接近於理想。要想辦到此事，矩陣的資料結構必須規劃成維度必須存在於資料結構中，求 1/[A] 時，此函數指令才能自動到資料結構中取得維度數據，然後才處理出結果。原本在 Win32Forth 中發展出來的 ABC Forth 系統並沒有這樣的規劃，2012 年以後，我轉往 Linux 作業系統的環境發展系統，才開始規劃出這樣的矩陣資料結構。這樣做，在想要印出矩陣內容時也會更為方便，只要能取得維度數據，印出矩陣的執行迴路指標就可以被自動設定，一個印出矩陣的指令就可以適用於印出任何維度的矩陣。
:
\ ABC FORTH demo code: Inverse Matrix
\ Author: Ching-Tang Tseng, Hamilton, New Zealand
\ Date: 2009-06-14 (renewed for c.l.f posting on 2014-12-18)
\ Contact: ilikeforth@gmail.com
\ Reference: F.R.Rukdeschel,"BASIC Scientific Subroutines Vol.1',
\ BYTE/McGRAW-HILL, 1980, p.95, Matrix inversion subroutine (MATINV)
\ Gauss-Jordan Elimination
\ Matrix A is input. Matrix C is ouptput
20 VALUE D \ Unique setting in this program for maximum matrix Dimension
D D MATRIX U
D D MATRIX V
D D MATRIX W
D D MATRIX A
D D 2 * MATRIX B
D D MATRIX C
7 INTEGERS I J K N I+N J+N M
1 REALS B1
((
: InputData ( -- )
[[ N = 5 ]]
{{ A ( 1 1 ) = 1 }} {{ A ( 1 2 ) = 0 }} {{ A ( 1 3 ) = 0 }} {{ A ( 1 4 ) = 0 }} {{ A ( 1 5 ) = 0 }}
{{ A ( 2 1 ) = 1 }} {{ A ( 2 2 ) = 1 }} {{ A ( 2 3 ) = 0 }} {{ A ( 2 4 ) = 0 }} {{ A ( 2 5 ) = 0 }}
{{ A ( 3 1 ) = 1 }} {{ A ( 3 2 ) = 1 }} {{ A ( 3 3 ) = 1 }} {{ A ( 3 4 ) = 0 }} {{ A ( 3 5 ) = 0 }}
{{ A ( 4 1 ) = 1 }} {{ A ( 4 2 ) = 1 }} {{ A ( 4 3 ) = 1 }} {{ A ( 4 4 ) = 1 }} {{ A ( 4 5 ) = 0 }}
{{ A ( 5 1 ) = 1 }} {{ A ( 5 2 ) = 1 }} {{ A ( 5 3 ) = 1 }} {{ A ( 5 4 ) = 1 }} {{ A ( 5 5 ) = 1 }}
;
))
: InputData ( -- )
[[ N = 5 ]]
{{ A ( 1 1 ) = 1 }} {{ A ( 1 2 ) = 3 }} {{ A ( 1 3 ) = 2 }} {{ A ( 1 4 ) = 6 }} {{ A ( 1 5 ) = 4 }}
{{ A ( 2 1 ) = 3 }} {{ A ( 2 2 ) = 5 }} {{ A ( 2 3 ) = 4 }} {{ A ( 2 4 ) = 2 }} {{ A ( 2 5 ) = 6 }}
{{ A ( 3 1 ) = 3 }} {{ A ( 3 2 ) = 8 }} {{ A ( 3 3 ) = 6 }} {{ A ( 3 4 ) = 7 }} {{ A ( 3 5 ) = 5 }}
{{ A ( 4 1 ) = 2 }} {{ A ( 4 2 ) = 5 }} {{ A ( 4 3 ) = 3 }} {{ A ( 4 4 ) = 4 }} {{ A ( 4 5 ) = 4 }}
{{ A ( 5 1 ) = 7 }} {{ A ( 5 2 ) = 6 }} {{ A ( 5 3 ) = 8 }} {{ A ( 5 4 ) = 7 }} {{ A ( 5 5 ) = 6 }}
;
: MatrixInverse ( -- )
BASIC
10 FOR I = 1 TO N
20 FOR J = 1 TO N
30 LET J+N = J + N
40 LET { B ( I J+N ) = 0 }
50 LET { B ( I J ) = A ( I J ) }
60 NEXT J
70 LET I+N = I + N
80 LET { B ( I I+N ) = 1 }
90 NEXT I
\ Perform row oeiented operations to convert the left hand side of B
\ to the identity matrix. The inverse of A will then be on the right.
100 FOR K = 1 TO N
110 IF K = N THEN 240
120 LET M = K
\ Find maximum element
130 FOR I = K + 1 TO N
140 IF { ABS ( B ( I K ) ) > ABS ( B ( M K ) ) } THEN 160
150 GOTO 170
160 LET M = I
170 NEXT I
180 IF M = K THEN 240
190 FOR J = K TO 2 * N
200 LET { B1 = B ( K J ) }
210 LET { B ( K J ) = B ( M J ) }
220 LET { B ( M J ) = B1 }
230 NEXT J
\ Divide row K
240 FOR J = K + 1 TO 2 * N
250 LET { B ( K J ) = B ( K J ) / B ( K K ) }
260 NEXT J
270 IF K = 1 THEN 340
280 FOR I = 1 TO K - 1
290 FOR J = K + 1 TO 2 * N
300 LET { B ( I J ) = B ( I J ) - B ( I K ) * B ( K J ) }
310 NEXT J
320 NEXT I
330 IF K = N THEN 350 \ Beware of this branch will be out of K loop.
340 GOTO 370
350 RUN 2drop
360 GOTO 430
370 FOR I = K + 1 TO N
380 FOR J = K + 1 TO 2 * N
390 LET { B ( I J ) = B ( I J ) - B ( I K ) * B ( K J ) }
400 NEXT J
410 NEXT I
420 NEXT K
\ Retrieve inverse from the right side of B
430 FOR I = 1 TO N
440 FOR J = 1 TO N
450 LET J+N = J + N
460 LET { B ( I J ) = B ( I J+N ) }
470 NEXT J
480 NEXT I
\ B is not a square matrix so transfer to square matrix C
490 FOR I = 1 TO N
500 FOR J = 1 TO N
510 LET { C ( I J ) = B ( I J ) }
520 NEXT J
530 NEXT I
490 END
;
: PrintMatrixW ( -- )
BASIC
10 RUN CR
20 FOR I = 1 TO N
30 FOR J = 1 TO N
40 LET { B1 = W ( I J ) } \ Based on W matrix
40 RUN B1 (F,) \ Right justified printing
50 NEXT J
60 RUN CR
70 NEXT I
80 END
;
: TotalAmount ( -- )
D 1+ DUP * 3 CELLS * ;
: MoveAtoW ( -- )
\ address-of A 0 0 FMatrix
\ address-of W 0 0 FMatrix
{{ A 0 0 }} {{ W 0 0 }} \ or
TotalAmount CMOVE ;
: MoveCtoW ( -- )
address-of C 0 0 FMatrix
address-of W 0 0 FMatrix
TotalAmount CMOVE ;
: MoveAtoU ( -- )
address-of A 0 0 FMatrix
address-of U 0 0 FMatrix
TotalAmount CMOVE ;
: MoveCtoV ( -- )
address-of C 0 0 FMatrix
address-of V 0 0 FMatrix
TotalAmount CMOVE ;
: Matrix* ( -- )
BASIC
\ Square Matrix Multiplication
\ W = U x V
10 FOR I = 1 TO N
20 FOR J = 1 TO N
30 LET { W ( I J ) = 0 }
40 FOR K = 1 TO N
50 LET { W ( I J ) = W ( I J ) + U ( I K ) * V ( K J ) }
60 NEXT K
70 NEXT J
80 NEXT I
90 END
;
: Checking ( -- )
BASIC
10 RUN MoveAtoU
20 RUN MoveCtoV
30 RUN Matrix*
40 RUN PrintMatrixW
50 END
;
: Once ( -- )
BASIC
10 PRINT " Matrix A = "
20 RUN MoveAtoW
30 RUN PrintMatrixW
40 RUN MatrixInverse
50 PRINT " Inverse of A = "
60 RUN MoveCtoW
70 RUN PrintMatrixW
80 PRINT " Matrix A times inverse of A = Idendity Matrix "
90 RUN Checking
100 END
;
: test ( -- )
BASIC
10 RUN InputData
20 RUN Once
30 END
;
: hi ( -- )
BASIC
10 PRINT " Input dimension of Matrix to be inverted "
20 INPUTI N
30 FOR I = 1 TO N
40 PRINT " Input row " ; I
50 FOR J = 1 TO N
60 RUN CR
70 INPUTR B1
80 LET { A ( I J ) = B1 }
90 NEXT J
100 NEXT I
110 RUN Once
120 END
;
: ReadData ( -- )
BASIC
10 RUN S" (38-1)MatInvData01.f" Get-File
20 RUN GetOneLineData
30 LET N = INT ( BRA ( 1 ) )
40 FOR I = 1 TO N
50 RUN GetOneLineData
60 FOR J = 1 TO N
70 LET { A ( I J ) = BRA ( J ) }
80 NEXT J
90 NEXT I
100 END
;
: main ( -- )
BASIC
10 RUN ReadData
20 RUN Once
30 END
;
page
.( Usage: ) cr
.( 1. test : for fixed input data ) cr
.( 2. hi : for interactive input data ) cr
.( 3. main : for data input from file ) cr
\S
Usage:
1. test : for fixed input data
2. hi : for interactive input data
ok
test
Matrix A =
1.00000 3.00000 2.00000 6.00000 4.00000
3.00000 5.00000 4.00000 2.00000 6.00000
3.00000 8.00000 6.00000 7.00000 5.00000
2.00000 5.00000 3.00000 4.00000 4.00000
7.00000 6.00000 8.00000 7.00000 6.00000
Inverse of A =
-.262431 -.367403 -.602210 1.17680 .259669
-.279006 -.190608 -.071823 .745856 -.060773
.135359 .389503 .668508 -1.48066 -.049724
.157459 -.179558 -.038674 .093923 .044199
.220994 .309392 -.071823 -.254144 -.060773
Matrix A times inverse of A = Idendity Matrix
1.00000 -2.22045E-16 -1.11022E-16 -2.22045E-16 -2.77556E-17
2.22045E-16 1.00000 -2.77556E-16 -2.22045E-16 5.55112E-17
0.00000E0 -2.22045E-16 1.00000 4.44089E-16 5.55112E-17
1.11022E-16 -2.22045E-16 -3.33067E-16 1.00000 -5.55112E-17
2.22045E-16 -8.88178E-16 -1.11022E-16 8.88178E-16 1.00000
ok
hi
Input dimension of Matrix to be inverted ? 3
Input row 1
? 1
? 0
? 0
Input row 2
? 0
? 2
? 0
Input row 3
? 0
? 0
? 3
Matrix A =
1.00000 0.00000E0 0.00000E0
0.00000E0 2.00000 0.00000E0
0.00000E0 0.00000E0 3.00000
Inverse of A =
1.00000 0.00000E0 0.00000E0
0.00000E0 .500000 0.00000E0
0.00000E0 0.00000E0 .333333
Matrix A times inverse of A = Idendity Matrix
1.00000 0.00000E0 0.00000E0
0.00000E0 1.00000 0.00000E0
0.00000E0 0.00000E0 1.00000
ok
3 sigdigits ! ok
main
Fileid is: 1212
139 Bytes are read into the Fadr RAM buffer!
Matrix A =
1.00 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0
1.00 1.00 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0
1.00 1.00 1.00 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0
1.00 1.00 1.00 1.00 0.00E0 0.00E0 0.00E0 0.00E0
1.00 1.00 1.00 1.00 1.00 0.00E0 0.00E0 0.00E0
1.00 1.00 1.00 1.00 1.00 1.00 0.00E0 0.00E0
1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.00E0
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
Inverse of A =
1.00 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0
-1.00 1.00 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0
0.00E0 -1.00 1.00 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0
0.00E0 0.00E0 -1.00 1.00 0.00E0 0.00E0 0.00E0 0.00E0
0.00E0 0.00E0 0.00E0 -1.00 1.00 0.00E0 0.00E0 0.00E0
0.00E0 0.00E0 0.00E0 0.00E0 -1.00 1.00 0.00E0 0.00E0
0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 -1.00 1.00 0.00E0
0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 -1.00 1.00
Matrix A times inverse of A = Idendity Matrix
1.00 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0
0.00E0 1.00 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0
0.00E0 0.00E0 1.00 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0
0.00E0 0.00E0 0.00E0 1.00 0.00E0 0.00E0 0.00E0 0.00E0
0.00E0 0.00E0 0.00E0 0.00E0 1.00 0.00E0 0.00E0 0.00E0
0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 1.00 0.00E0 0.00E0
0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 1.00 0.00E0
0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0 1.00
ok
The format of data to be stored in MatInvData01.f is as folloing:
8
1 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0
1 1 1 0 0 0 0 0
1 1 1 1 0 0 0 0
1 1 1 1 1 0 0 0
1 1 1 1 1 1 0 0
1 1 1 1 1 1 1 0
1 1 1 1 1 1 1 1
第(39)個範例是一個實用程式，我用它來計算出機械設計時購買傳動皮帶所需要的皮帶長度數據。
程式內使用幾何公式進行精確計算，購買時，就買最大長度略小於理論計算值的標準規格現貨皮帶。這個程式用起來屢試不爽，後來，我便專為隔壁不會用電腦的老鄰居 Mr. Norman 提供服務，他是擁有車床技術的老技工，採用繞繩量長度的經驗去買皮帶時，失誤過兩次，曾經拆用新品後退換，店家賣的有點不爽，他的經驗作法不如我的理論計算來得精準，一次就能搞定，於是，每次他要再買皮帶前，就找我用這個電腦程式計算出結果。
1996 年初，我剛到紐西蘭時，身上沒帶太多的錢來移民，沒錢買新機器。自己又很喜歡動手修理或製作工具，所以經常去二手店搜尋廉價的舊品，回家自己組裝出好用的工作母機。我熟各方面的科技，是個有合格證照的工程師，主修機械、懂電機與電子、敢搞化學與生物、熟悉許多的儀器、理論與實務兼備，更擅長設計電腦系統。因此，車庫內庫存的機器，都是我自己組裝出來的，當然，需要動用車床與電焊時，我只好找隔壁的老先生幫忙，他沒拒絕過，還特別喜歡定期找我搞新產品，我車庫內的巨型工作桌檯，就是他去專科學校做專業 Sharpen technician 時，學校淘汰下來轉送我的贈品，它厚重結實到可以使用一輩子。
這個範例，順便介紹涉及這個範例程式的圖解說明與產品實物，貼文附上真實照片，組裝這些機器時，都用過這個範例程式。
大學的機械設計課程中，皮帶是專門的一章，機械常要用到。這個範例程式，就算寫足了參考文獻資料來源，大家可能還是很難弄懂所用公式的對應意義，我用兩張自己準備的示意圖，來表示公式的理論依據，一看就能明白，一張是平行安裝的皮帶，另一張是交叉安裝的皮帶。
簡圖兩張
工作母機分成車、銑、鉋、磨、鑽、研、搪，我沒那麼多需求，最需要的則是鑽床。
照片一是這台鑽床的照片，就是用我廉價買回來的機架回收品組裝而成的，機架還算完整。
夾頭是個台灣朋友在桃園平鎮開鈕扣製造機械工廠的老闆所贈送的，鑽床原來的舊馬達燒掉了，皮帶也斷掉沒有了。
我花 5 元紐幣，買了個舊的半馬力感應馬達，畫了張草圖，麻煩隔壁老先生幫我車製了一個鋁合金的皮帶盤，自製了一片平板固定鉄板，鎖上螺絲固定馬達後，用卡尺量了上述貼圖中所需要的相關參數，就去買了這條新皮帶 19 元紐幣。
待試車運轉時，才發現轉動方向相反了。 Norman 老先生不懂電，我懂，我就拆開馬達的接線盒，量出啟動線圈與運轉線圈的接頭，做了標記，在書桌上畫出線路圖，確定反裝一組接頭就能令轉向倒過來，於是自已動手換了接線頭的位置。很險，四條中有一條幾乎拉到了極限才能鎖定，重新再試，成功了。於是把修正後的接線圖記進筆記本，永久留參，又從舊瓶罐中找出一個庫存能耐較大電流的扳動開關，裝機固定，就完成了。
這台鑽床，我已用了 30 幾年，常常用。後來回台灣，還去買了一個工件夾具，這樣，鑽床就也能銑出木料，產出簡單有花樣的東西了。
標準規格皮帶現貨，程式算完後一次就買定了，買時，我很有信心，買到就走，直接安裝使用，一直用到今天。照片中皮帶盒的蓋子是為了照相方便才掀開的。
照片一、鑽床
第二張照片，比較簡單，是鑽頭與砂輪同軸雙用的廢機器，廢品回收廠論斤秤重賣我 3 元紐幣，馬達還是好的，皮帶早就斷掉沒有了，砂輪軸長期沒有潤滑，早就咬死不能轉動，我運回家後，先加油潤滑，讓軸承能夠順利運轉，然後用電腦算出皮帶長，由於兩邊的滑輪尺寸一樣，皮帶比較好找，我就從廢車場的 yard 內拆買回了一條 1 元汽車引擎內發電機用的舊皮帶，長度比程式算出的結果稍長，只是相當，且日貨的規格與英制的規格有點差別，只好在安裝皮帶前，調整馬達的固定螺絲，調到極限了，還是有點不夠延伸，只好將就著、鬆鬆的用了。
我只用這台砂輪機來磨除金屬會傷手的毛邊，都只是粗磨，不做精確研磨使用，馬達原就是感應馬達，大概使用一輩子也燒不壞，加潤滑油的油孔蓋，有一個已經丟失了，只要常用、常加油，沒蓋子也無妨，我磨很髒的東西時，才會用一個小木塞塞住潤滑油加油孔。
照片二、砂輪機
第三張照片，是一台改裝後的桌台式圓鋸，由於精度全無，我大約一年才會用一次，用來鋸短庭院中的大樹幹，送進爐膛，冬天取暖時燒掉。
機器原本很好、也很貴，是本市一位開畫廊的台灣人老闆，誤用此機，鋸室內裝潢硬質材料時，過載使用，燒掉了並激式有碳刷馬達的碳刷，不宜換修了，才丟給我的。
我當然不會、也不可能再買一個新的同款很貴的馬達回裝了，但我換算出力需求與需要轉速後，算出一馬力的感應馬達能夠驅動。於是，去二手店花 15 元買了一台一馬力、單相、兩極、能快轉的大型舊的感應馬達，因為沒有碳刷才能再用一輩子。馬達當然不可能裝回原來空間很小的有碳刷馬達位置了，只好自行設計安裝支架、算出能讓轉速高達每分鐘 8000 轉左右的大小皮帶盤，一個大的買 1 元舊品，另一個小的請 Norman 老先生幫忙用車床產製新品，就裝成了。皮帶也是用電腦程式算出來的，一次買定，用到現在。
紐西蘭的電力是 240 伏特 50 Hz 的家用供電，理論上兩極的感應馬達，每分鐘最高轉速不會大於 3000 轉，一定會有相位落後的轉差，無負荷運轉也不會高過 3000 轉，皮帶盤變速三倍後約為 9000 轉，傳動磨擦損失後大約還能有 8000 轉，馬達原裝標示牌上說有一萬轉，我重裝，就轉慢一點， 9000 轉也沒有關係，若設計成轉得太快，怕會崩盤，輸出扭力夠，就行。
照片三、桌鋸的皮帶盤
:
程式內使用幾何公式進行精確計算，購買時，就買最大長度略小於理論計算值的標準規格現貨皮帶。這個程式用起來屢試不爽，後來，我便專為隔壁不會用電腦的老鄰居 Mr. Norman 提供服務，他是擁有車床技術的老技工，採用繞繩量長度的經驗去買皮帶時，失誤過兩次，曾經拆用新品後退換，店家賣的有點不爽，他的經驗作法不如我的理論計算來得精準，一次就能搞定，於是，每次他要再買皮帶前，就找我用這個電腦程式計算出結果。
1996 年初，我剛到紐西蘭時，身上沒帶太多的錢來移民，沒錢買新機器。自己又很喜歡動手修理或製作工具，所以經常去二手店搜尋廉價的舊品，回家自己組裝出好用的工作母機。我熟各方面的科技，是個有合格證照的工程師，主修機械、懂電機與電子、敢搞化學與生物、熟悉許多的儀器、理論與實務兼備，更擅長設計電腦系統。因此，車庫內庫存的機器，都是我自己組裝出來的，當然，需要動用車床與電焊時，我只好找隔壁的老先生幫忙，他沒拒絕過，還特別喜歡定期找我搞新產品，我車庫內的巨型工作桌檯，就是他去專科學校做專業 Sharpen technician 時，學校淘汰下來轉送我的贈品，它厚重結實到可以使用一輩子。
這個範例，順便介紹涉及這個範例程式的圖解說明與產品實物，貼文附上真實照片，組裝這些機器時，都用過這個範例程式。
大學的機械設計課程中，皮帶是專門的一章，機械常要用到。這個範例程式，就算寫足了參考文獻資料來源，大家可能還是很難弄懂所用公式的對應意義，我用兩張自己準備的示意圖，來表示公式的理論依據，一看就能明白，一張是平行安裝的皮帶，另一張是交叉安裝的皮帶。
工作母機分成車、銑、鉋、磨、鑽、研、搪，我沒那麼多需求，最需要的則是鑽床。
照片一是這台鑽床的照片，就是用我廉價買回來的機架回收品組裝而成的，機架還算完整。
夾頭是個台灣朋友在桃園平鎮開鈕扣製造機械工廠的老闆所贈送的，鑽床原來的舊馬達燒掉了，皮帶也斷掉沒有了。
我花 5 元紐幣，買了個舊的半馬力感應馬達，畫了張草圖，麻煩隔壁老先生幫我車製了一個鋁合金的皮帶盤，自製了一片平板固定鉄板，鎖上螺絲固定馬達後，用卡尺量了上述貼圖中所需要的相關參數，就去買了這條新皮帶 19 元紐幣。
待試車運轉時，才發現轉動方向相反了。 Norman 老先生不懂電，我懂，我就拆開馬達的接線盒，量出啟動線圈與運轉線圈的接頭，做了標記，在書桌上畫出線路圖，確定反裝一組接頭就能令轉向倒過來，於是自已動手換了接線頭的位置。很險，四條中有一條幾乎拉到了極限才能鎖定，重新再試，成功了。於是把修正後的接線圖記進筆記本，永久留參，又從舊瓶罐中找出一個庫存能耐較大電流的扳動開關，裝機固定，就完成了。
這台鑽床，我已用了 30 幾年，常常用。後來回台灣，還去買了一個工件夾具，這樣，鑽床就也能銑出木料，產出簡單有花樣的東西了。
標準規格皮帶現貨，程式算完後一次就買定了，買時，我很有信心，買到就走，直接安裝使用，一直用到今天。照片中皮帶盒的蓋子是為了照相方便才掀開的。
第二張照片，比較簡單，是鑽頭與砂輪同軸雙用的廢機器，廢品回收廠論斤秤重賣我 3 元紐幣，馬達還是好的，皮帶早就斷掉沒有了，砂輪軸長期沒有潤滑，早就咬死不能轉動，我運回家後，先加油潤滑，讓軸承能夠順利運轉，然後用電腦算出皮帶長，由於兩邊的滑輪尺寸一樣，皮帶比較好找，我就從廢車場的 yard 內拆買回了一條 1 元汽車引擎內發電機用的舊皮帶，長度比程式算出的結果稍長，只是相當，且日貨的規格與英制的規格有點差別，只好在安裝皮帶前，調整馬達的固定螺絲，調到極限了，還是有點不夠延伸，只好將就著、鬆鬆的用了。
我只用這台砂輪機來磨除金屬會傷手的毛邊，都只是粗磨，不做精確研磨使用，馬達原就是感應馬達，大概使用一輩子也燒不壞，加潤滑油的油孔蓋，有一個已經丟失了，只要常用、常加油，沒蓋子也無妨，我磨很髒的東西時，才會用一個小木塞塞住潤滑油加油孔。
第三張照片，是一台改裝後的桌台式圓鋸，由於精度全無，我大約一年才會用一次，用來鋸短庭院中的大樹幹，送進爐膛，冬天取暖時燒掉。
機器原本很好、也很貴，是本市一位開畫廊的台灣人老闆，誤用此機，鋸室內裝潢硬質材料時，過載使用，燒掉了並激式有碳刷馬達的碳刷，不宜換修了，才丟給我的。
我當然不會、也不可能再買一個新的同款很貴的馬達回裝了，但我換算出力需求與需要轉速後，算出一馬力的感應馬達能夠驅動。於是，去二手店花 15 元買了一台一馬力、單相、兩極、能快轉的大型舊的感應馬達，因為沒有碳刷才能再用一輩子。馬達當然不可能裝回原來空間很小的有碳刷馬達位置了，只好自行設計安裝支架、算出能讓轉速高達每分鐘 8000 轉左右的大小皮帶盤，一個大的買 1 元舊品，另一個小的請 Norman 老先生幫忙用車床產製新品，就裝成了。皮帶也是用電腦程式算出來的，一次買定，用到現在。
紐西蘭的電力是 240 伏特 50 Hz 的家用供電，理論上兩極的感應馬達，每分鐘最高轉速不會大於 3000 轉，一定會有相位落後的轉差，無負荷運轉也不會高過 3000 轉，皮帶盤變速三倍後約為 9000 轉，傳動磨擦損失後大約還能有 8000 轉，馬達原裝標示牌上說有一萬轉，我重裝，就轉慢一點， 9000 轉也沒有關係，若設計成轉得太快，怕會崩盤，輸出扭力夠，就行。
:
\ PulleyBeltLength.f
((
Theme : The length of pulley's belt
Author : Ching-Tang Tseng, Hamilton, New Zealand
Date : 2015-01-18
Reference : Frank Ayers, JR."Plane and spherical Trigonometry",
Schaum publishing Co., 1954, p.115
Description :
Find the length (1)L1 and (2)L2 of
(1)an open driving belt
(2)a crossed driving belt
running around two pulleys of radii B, and s, respectively
if the distance between the centres of the wheels is d
The formula are:
Theta T1 = Arcsin((B-s)/d), then
Length L1 = (B+s)*pi + 2*(B-s)*T1 + 2*d*cos(T1)
e.g. B=15, s=5, d=30 then L1=126.2
Theta T2 = Arcsin((B+s)/d), then
Length L2 = (B+s)*(pi+2*T2) + 2*(B-s)*T2 + 2*d*cos(T2)
e.g. B=10, s=5, d=30 then L2=114.8
))
7 Reals B s d T1 L1 T2 L2
: hi ( -- )
basic
\ QAEnter data
10 print " Enter Big R, Small r, Distance d ( B s d )="
20 inputr B , s , d
\ Enter data rational checking
30 if { d >= ( B + s ) } then 60
40 print " Warning! d must >= (B+s), Re-Enter Please."
50 goto -10
\ Calculation by formula
60 let { T1 = ASIN ( ( B - s ) / d ) }
70 let { L1 = ( B + s ) * fpi + 2 * ( B - s ) * T1 + 2 * d * cos ( T1 ) }
80 let { T2 = ASIN ( ( B + s ) / d ) }
90 let { L2 = ( B + s ) * ( fpi + 2 * T2 ) + 2 * d * cos ( T2 ) }
\ Results output
100 print " The length of an open driving belt is : " ; { L1 }
110 print " The length of a crossed driving belt is : " ; { L2 }
120 end ;
\S
hi
Enter Big R, Small r, Distance d ( B s d )=? 15 5 30
The length of an open driving belt is : 126.197
The length of a crossed driving belt is : 136.742 ok
hi
Enter Big R, Small r, Distance d ( B s d )=? 10 5 30
The length of an open driving belt is : 107.959
The length of a crossed driving belt is : 114.793 ok
第(40)個範例，是一個很簡單的線性插值程式。
數值分析教科書的第一章，通常講的就是插值的數值分析算法，方法很多，技術也一直隨著時代的進步而發展。這一個範例是線性的求法，工程應用上也很有用，教科書上反而不講，教學都集中在非線性的各種插值法的理論探討上，包括內插、外插、勻插、不規則分割距離的插值運用.....等等。
在工程實務上，插值才是最逼真的數據，技術雖很基礎卻很有用，人類最複雜的電腦應用例如：分析氣象資料。氣象資料的收集都不可能勻稱，於是就得依搞大量的插值運算進行數據處理，後續複雜的分析，才能再根據全部的插值結果算出梯度，然後才能預測天氣、風向、晴雨、繪製動態氣象預測圖...等等，做這些事情，都在搞大量的插值計算。
我寫的範例程式很簡單，而且採用了中文的對答式輸入法進行計算，也做了不合理數據輸入時的防止設計，請載入程式後執行之，就能體會出結果。如果我用純粹的 Forth 格式寫程式，保證大家一定看不懂在做什麼，這就是此範例簡單而特別之處。
多維、不等距離、定範圍內的插值技術採用差分法求解答案，基本原則就是把間隔分得越細，所得結果的趨勢變化就能越平滑，但計算就越耗用記憶體也耗用計算時間，近代電腦已不在乎這些消耗。20150413 貼出的網文： Curve fitting demo in ABC FORTH 有這樣的例子，程式中宣告了 57 層的陣列供差分表使用，演算結果繪製成曲線，顯示出它的逼近度，也就是表示效果已經夠好到了什麼程度。
:
數值分析教科書的第一章，通常講的就是插值的數值分析算法，方法很多，技術也一直隨著時代的進步而發展。這一個範例是線性的求法，工程應用上也很有用，教科書上反而不講，教學都集中在非線性的各種插值法的理論探討上，包括內插、外插、勻插、不規則分割距離的插值運用.....等等。
在工程實務上，插值才是最逼真的數據，技術雖很基礎卻很有用，人類最複雜的電腦應用例如：分析氣象資料。氣象資料的收集都不可能勻稱，於是就得依搞大量的插值運算進行數據處理，後續複雜的分析，才能再根據全部的插值結果算出梯度，然後才能預測天氣、風向、晴雨、繪製動態氣象預測圖...等等，做這些事情，都在搞大量的插值計算。
我寫的範例程式很簡單，而且採用了中文的對答式輸入法進行計算，也做了不合理數據輸入時的防止設計，請載入程式後執行之，就能體會出結果。如果我用純粹的 Forth 格式寫程式，保證大家一定看不懂在做什麼，這就是此範例簡單而特別之處。
多維、不等距離、定範圍內的插值技術採用差分法求解答案，基本原則就是把間隔分得越細，所得結果的趨勢變化就能越平滑，但計算就越耗用記憶體也耗用計算時間，近代電腦已不在乎這些消耗。20150413 貼出的網文： Curve fitting demo in ABC FORTH 有這樣的例子，程式中宣告了 57 層的陣列供差分表使用，演算結果繪製成曲線，顯示出它的逼近度，也就是表示效果已經夠好到了什麼程度。
:
\ (40)LinearInterpolation.f
((
Theme : Linear interpolation
Author : Ching-Tang Tseng, Hamilton, New Zealand
Date : 2015-01-25
Reference : NAN
Description :
Yi = Y1+(Xi-X1)*(Y2-Y1)/(X2-X1), X1<X2
))
6 Reals X1 Y1 X2 Y2 Xi Yi
: RESET
{{ X1 = 0 }} {{ Y1 = 0 }}
{{ X2 = 0 }} {{ Y2 = 0 }}
{{ Xi = 0 }} {{ Yi = 0 }}
;
: QAEnter BASIC
10 RUN PAGE CR
20 PRINT " 線性內插程式，根據兩組已知之座標值，計算出位於其間的內插值。 "
30 PRINT " 本程式要求輸入之兩點座標之值，X1 必須小於 X2，請注意 X1<X2！ "
40 RUN CR CR
50 PRINT " 請輸入第一組已知數據 X1 , Y1 : "
60 InputR X1 , Y1
70 PRINT " 請輸入第二組已知數據 X2 , Y2 : "
80 InputR X2 , Y2
90 END
;
Integer K
: hi BASIC
10 RUN RESET
20 RUN QAEnter
30 PRINT { " 待計算插值之 Xi 座標，被限制為介於 " ; X1 ; " 與 " ; X2 ; " 之間 " }
40 PRINT " 請輸入 Xi 之值 "
50 InputR Xi
60 IF { ( Xi >= X1 ) AND ( Xi <= X2 ) } THEN 90
70 PRINT " Xi 超出範圍，請重新輸入！ "
80 GOTO -30
90 LET { Yi = Y1 + ( Xi - X1 ) * ( Y2 - Y1 ) / ( X2 - X1 ) }
100 PRINT { " 當 Xi = " ; Xi ; " 時 " ; " Yi = " ; Yi }
110 PRINT CR " 繼續計算 Xi 之 Yi 值，Repeat(Y/N)? "
120 LET K = KEY
130 IF ( K = 89 ) OR ( K = 121 ) THEN -40
140 END
;
訂閱：
張貼留言 (Atom)
沒有留言:
張貼留言