前言:
是因為這樣的, 在工作上, 台灣過去20年來培養了大量的製程工程師. 有製程工程師就會有大量的工程實驗, 或是製程實驗. 有製程實驗就免不了產生大量的實驗數據要分析, 要歸納. 於是, 大家都一定不陌生的是最小平方法的曲線配置. 然後, 每個工程師們都不陌生就是如何用 Excel 來做這個最小平方法的多項式配置.
這邊, 不甘於只是用 Excel 來做這件事啦! 想說自己來寫一個程式來做這個最小平方法的多項式配適 (這個程式已經完成了, 跟這篇 BLOG 剛好是姊妹篇, 有興趣請參照此連結), 這樣一定很方便. 如果你熟悉數值分析技術的話會知道, 多階的最小平方法多項式曲線配置, 最後會轉化變成要解一個n階的線性聯立方程組問題. 所以這邊第一步要把這個可以解任意階的聯立方程組的數值程式給先搞定. 所以就來寫一個吧, 來寫一個可以解 n 階聯立方程組的程式囉!
我覺得啊, 數值分析也是工程師們必備的技能. 而且, 結合數學, 用電腦演算法來計算, 模擬工程或物理問題, 或用數理統計處理實驗數據, 這一直都是非常有趣的事. 當然啦, 為了推廣我所非常喜愛的很稀有的 FORTH 程式語言, 這裡的程式一概都用 FORTH 語言來撰寫(C語言, 或其他語言的數值分析程式, 真的太普遍了, 不差我這一個.). 順便展現一下 FORTH 語言優秀的高階性, 跟交談性, 及用來執行數值分析的技巧. 當然, 這些程式碼也變成我自己以後常用的數值分析程式庫啦!
FORTH 科學計算:
講到科學計算, 一定離不開浮點運算. FORTH語言的發明人, 祖師爺 Chuck Moore 其實本身蠻排斥浮點運算的! 因為早期的電腦蠻不威的, 大部分的 CPU 都沒有浮點運算的硬體. 強用浮點運算的結果就是雖然方便, 但是速度會變得很慢. 所以死硬派的 FORTH 程式設計師是不用浮點運算的.
神級般的祖師爺 Chuck Moore, 也身體力行, 他用他發明的新一代 FORTH 語言, color FORTH, 在完全不用浮點運算的情況下, 只用固定整數運算就開發出 OKCAD 專門在做晶片設計開發模擬的神級的 CAD 軟體. 真是難以置信這類高度依賴物理定律來計算裡頭數萬數億顆電晶體電壓電流及溫度運作情況的 CAD 軟體可以不依賴浮點運算來計算跟運作.
Starting FORTH 其實也做了很好的展示, 很多情況我們認為非浮點運算不可的計算, 其實固定整數運算也可以作, 而且做得更好, 速度更快. 所以歷屆主流的 FORTH 程式設計師是不太用符點運算的. 所以早期的 FORTH, 這塊爭議最大, 大家用的方法也一團亂.
但是, 對於科學工作者, 浮點運算還是有很大的方便性. 畢盡, 我們不是聰明絕頂如神般的 FORTH 祖師爺 Chuck Moore, 可以自由駕馭整數運算或是精心打造如神般的軟體. 我們都是凡人, 浮點運算還是有一定的必要性.
於是直到 ANSI FORTH 美國國家標準, 對於 FORTH 浮點運算這塊, 開始做了一些標準的制定. 這樣的亂象才開始逐漸統一.
我們不是祖師爺, 對於現在的浮點運算能力超強的微處理器, 浮點運算不用白不用. 而且 FORTH 本身的高階性, 配上 ANSI FORTH 已經把這塊標準制定下來了, 這使得 FORTH 變得一個非常適合做數值分析處理的強大工具. 科學函式庫與其使用那些最後看起來像鬼畫符, 無法單獨交談測試的 FORTRAN或C, 我覺得可以單獨交談測試, 可以任意打造非常高階的語法, 不管高階低階都游刃有餘的 FORTH 會是更佳的選擇.
這樣的邏輯下, 已經有 FORTH 先驅者幫 FORTH 開發了一系列的 Open Source 的科學程式函式庫, 也被眾多的 FORTH 版本採用. 有興趣的大家可以參考喲. 寫得很巧妙. 但仔細研究過程式碼後, 我其實不是很喜歡這一系列程式庫的寫法啦, 一點 FORTH 味道都沒有, 搞得像是 C 和 FORTRAN 的綜合體.
還是自己來寫跟打造自己所要的比較對味些!
陣列的資料結構:
這裡要做的是利用高斯消去法來解聯立方程組, 這個數值方法是以數學陣列運算來當基礎的. 然後, 之後很多的數值方法也都是以數學陣列運算為基礎的. 所以為了兼顧之後使用的彈性, 最快的方式還是實作一個高階陣列吧!
先來設計它的語法:
5 7 DefMatrix Test 定義一個有 5列 (Row) 7欄 (Column) 名字叫做 Test 的陣列. Row=0..4, Column=0..6
Test RowNum@ 取出 Test 陣列裡列 (Row) 的總數, 這裡會傳回 5
Test ColNum@ 取出 Test 陣列裡欄 (Column) 的總數, 這裡會傳回 7
3 4 Test Matrix f@ 取出 Test 陣列裡(Row=3, Col=4)的內容(浮點數)
3.14e0 0 2 Test Matrix f! 將 3.14 這個浮點數存入 Test 陣列裡(Row=0, Col=2)的位置.
Test EmptyMatrix 將 Test 陣列裡的內容全部清除. (填零)
Test .matrix 列印 Test 陣列裡的內容
資料結構設計很簡單的啦, 前兩個 cell 存 欄(Column) 跟 列(Row) 的數目, 後面存浮點數的陣列內容. 完全手動創造的結構就像下面這樣.
create test1
4 5 , ,
1.1e0 f, 7.2e0 f, 1.3e0 f, 1.4e0 f, 1.5e0 f,
2.1e0 f, 2.2e0 f, 7.3e0 f, 2.4e0 f, 4.5e0 f,
1.1e0 f, -3.2e0 f, 7.3e0 f, 3.5e0 f, -2.0e0 f,
4.1e0 f, 4.5e0 f, 5.1e0 f, -4.4e0 f, 7.5e0 f,
所以, DefMatrix
: DefMatrix
create ( RowNum ColNum --)
2dup , , * floats allot
;
用 create 造個詞的頭, 然後把 欄(Column) 跟 列(Row) 兩個數字用 ',' 編進字典, 然後用 allot 為這個詞預留 RowNum*ColNum 個 floats 的浮點數空間下來.
資料結構的第一個 cell 存放的是陣列裡的 總欄(Column) 數目, 第二個 cell 存放的是陣列裡的 總列(Row) 數目,
所以 RowNum@ 跟 ColNum@
: ColNum@ ( maddr -- ColNum)
@
;
: RowNum@ ( maddr -- RowNum)
cell+ @
;
陣列內 第 (Row, Col) 位置的真實儲存的記憶位址的計算, 公式 addr = maddr + 2 cells + (Col + ColNum*Row) floats
用了一些堆疊運算來實作, 常用的 FORTH 技巧, 把擋路的數字利用 >r 推到返回堆疊暫存, 利用 r@ 複製回來使用, 非常方便. 最後記得 r> 再從返回堆疊推回參數堆疊, 一切船過了無痕!
: Matrix ( Row Col maddr -- addr')
>r r@ @ ( Row Col ColNum R: maddr)
rot * ( Col ColNum*Row R: maddr)
+ floats [ 2 cells ] literal +
r> +
;
當然啦, FORTH 最威的地方在這裡, 日後, 假如發現這個地方是程式執行速度的瓶頸. 可以把 FORTH 對這個微處理器的組合語言編譯器詞彙叫出來, 直接線上編譯把這段低階化, 取得最大執行速度的改善.
Dummy Matrix 指標變數,
這裡為了方便, 我們利用 value 來定義變數, 用來指向我們正要操作的陣列位址. (C語言稱這個叫做指標, FORTH 裡不管這樣的稱呼的, 只是個記憶體的位址而已沒必要弄得神秘兮兮的), 所以
: DummyMatrix
0 value
;
FORTH 裡有兩種定義變數的方法, 一個是 variable, 另一個是 value. 這裡採用 value, 其定義出來的變數是用 to 這個指令來改變裡頭的值的, 所以操作方式如下喲
DummyMatrix MatSelected 定義一個叫做 MatSelected 的變數 (這裡的用途叫做指標變數)
4 5 DefMatrix ExperimentData 定義一個 4x5大小 叫做 ExperimentData 的陣列
ExperimentData to MatSelected 將 ExperimentData 的陣列位址, 存入名為 MatSelected 的變數 (指標變數)
所以現在 MatSelected 指向 ExperimentData, 所以下面這兩個操作的結果會一樣.
ExperimentData .matrix 列印 ExperimentData 陣列內容
MatSelected .matrix 列印 MatSelected 所指向的陣列的陣列內容 (這裡已經指向 ExperimentData)
所以只要利用 to, 可以很方便的改變 MatSlelectd 這個 dummyMatrix 所定義出來的指標變數的內容, 指向不同想要使用的陣列.
來開始解說高斯消去法解聯立方程組囉!
三個個重要的數學規則如下:
1. 陣列裡列(Row)跟列(Row)的交換, 不會影響解的位置或數值.
2. 某一列整列的元素乘上一個常數, 不會影響解的位置或數值.
3. 某一列整列的元素乘上一個常數, 加上另外一列, 不會影響解的位置或數值.
解說:
1. 陣列裡列(Row)跟列(Row)的交換, 不會影響解的位置或數值.
很簡單呀, 舉例:
2 x1 + x2 - x3 = 8
-3 x1 - x2 + 2 x3 = -11
-2 x1 + x2 + 2 x3 = -3
對調前面兩列,
-3 x1 - x2 + 2 x3 = -11
2 x1 + x2 - x3 = 8
-2 x1 + x2 + 2 x3 = -3
所有的式子都沒變, 所以什麼都不會變的. XD
程式碼
給 Row1, Row2 要對調的兩行的參數, 用個 do-loop 把這兩行的每一個欄元素對調囉!
: RowSwap ( row2 row1 --)
2dup = if 2drop exit then 兩列要對調是同一行Row1=Row2,同一行沒什麼好調的,所以退出.
MatSelected ColNum@ 0 掃過每一欄 Col 元素的 do-loop
do over i MatSelected Matrix dup f@ 將row2的欄元素col=i的陣列值取出,順便留下位址
over i MatSelected Matrix dup f@ 將row1的欄元素col=i的陣列值取出,順便留下位址
fswap f! f! 交換一下內容, 利用剛剛留下的位址, 回存回陣列
loop 2drop
;
2. 某一列整列的元素乘上一個常數, 不會影響解的位置或數值.
一般代數的規則, 等號兩邊乘上一個常數, 兩邊還是相等地! , 舉例:
-3 x1 - x2 + 2 x3 = -11
給他們乘上 10, 兩邊還是相等地
(-3 x1 - x2 + 2 x3 = -11) * 10 相等 -30 x1 - 10 x1 + 20 x3 = -110
程式碼
遵循 FORTH 後序的數字計算習慣, Row*Scalar 應該寫成 Row Scalar * 所以將矩陣中某列呈上一個常數 Scalar 的指令命名為 RowScalar*
也是很簡單, 用個 do-loop 掃過指定 row 的那列的每個陣列元素, 取出乘上常數 scalar 後回存
: RowScalar* ( row F: scalar --)
MatSelected ColNum@ 0
do dup i MatSelected Matrix dup f@ fover f* f!
loop drop
;
3. 某一列整列的元素乘上一個常數, 加上另外一列, 不會影響解的位置或數值.
這還是來自於一般代數的規則,
2 x1 + x2 - x3 = 8 ----(1)
-3 x1 - x2 + 2 x3 = -11 ----(2)
-2 x1 + x2 + 2 x3 = -3
將 (1) * 3 + (2) -> (2) 所有的等式還是不會變呀! (等號兩邊同時加上一樣的東西, 結果還是一樣)
2 x1 + x2 - x3 = 8
3 x1 + 2 x2 - x3 = 13
-2 x1 + x2 + 2 x3 = -3
程式碼
一樣, 根據 FORTH 的運算規則, 取名為 Row2 Row1 Scalar * - to Row2; Row2Row1Scalar*Subtract>Row2, 給定 row2, row1, scalar 的數值, 會將 Row1 的所有陣列元素乘上 scalar 後跟 Rwo2 相減後存回 row2
也是不困難, 用個 do-loop 掃過指定 row2, row1 兩列的每個陣列元素, row1乘上常數 scalar 後和 row1相減後回存
: Row2Row1Scalar*Subtract>Row2 ( Row2 Row1 F: scalar --)
MatSelected ColNum@ 0
do over i MatSelected Matrix dup f@ ( r2 r1 r2adr F: sc r2)
fover ( r2 r1 r2adr F: sc r2 sc)
over i MatSelected Matrix f@ f* ( r2 r1 r2adr F: sc r2 r1*sc)
f- f!
loop 2drop fdrop
;
高斯消去法解聯立方程組
整個演算法規則如下:
聯立方程組,
等於這樣的矩陣,
只要想辦法將這個矩陣變成這樣, 只剩下中間斜對角這排, 且都是1 (單位矩陣), 自然就解到答案了, 因為此時 x1 = r''1, x2 = r''2, x3 = r''3, x4 = r''4
要怎麼將矩陣變成這樣呢?
先將下半部三角形的地方全部想辦法消去變成零呀, 然後再將上半部三角形的地方全部想辦法消去變成零呀, 這樣就是囉!
所以,
: SolveLinearEQ ( maddr --)
to MatSelected 將陣列位址放到指標變數,方便所有的指令來操作
AdjustPivot 每一欄尋找最大絕對值來當主軸,以減少浮點運算的 round-off 誤差
LowerTrianglize 利用每列主軸元素,執行下三角消去
UpperTrianglize 利用每列主軸元素,執行上三角消去
;
如何將下半部三角形的地方消去變成零呢?? LowerTrianglize
: LowerTrianglize
MatSelected RowNum@ 0
do
MatSelected RowNum@ i
do
i j MatSelected Matrix f@ f0<>
陣列元素已經為零的話, 什麼都不用做!
if
i j RowUnitizeByCol
將 i Row整列,以 j Col的數字來除,這樣可以讓 i Row的 j Col(主軸)的元素變成一
i j =
if j j MatSelected Matrix f@ f0= abort" Pivot is ZERO!!"
主軸為零, 無解!!! Abort吧!
else i j 1e0 Row2Row1Scalar*Subtract>Row2 then
將 i Row 整列跟 j Row 整列相減後放回 i Row
then
loop
loop
;
兩個 do-loop, 外面那個 掃過每個主軸欄位 Column, 內部內個 掃過每個主軸的下方列 Row, 形成掃過下面的三角形的所有元素.
指令 RowUnitizeByCol 先將第 i 列的所有元素, 以 j 欄的數字來除, 這樣可以讓第j欄(主軸)的數字全部都變成1. (根據規則2. 某一列整列的元素乘上一個常數 Scalar, 不會影響解的位置或數值.)
再來利用已經變成1 的 j Row 主軸列, 來消去其他主軸下方的所有 i Row 非主軸列. (根據規則3. 某一列整列的元素乘上一個常數, 加上另外一列, 不會影響解的位置或數值.)
這裡有兩個作法:
方法一: (非主軸列先除變成一後再減主軸列相減消去)
利用 RowUnitizeByCol 將 j Row主軸及其下方先變成 1, 再用 Row2Row1Scalar*Subtract>Row2 將它們相減, 這樣 j Row 主軸下方就消去變成零囉!
所以, i RowUnitizeBy j Col 後, 再用 i Row2 j Row1 1.0e1 f* 放入 i Row2 (指令 i j 1e0 Row2Row1Scalar*Substract>Row2)
方法二: (主軸列為一的值直接乘上非主軸列相等後再相減消去)
利用 RowUnitizeByCol 將 i Row主軸變成 1 (需放在外迴圈, 指令 i i RowUnitizeByCol)
內迴圈內, 再將主軸列 j Row 乘上 i Row 於主軸 j Col 的數字值, 相減就減掉囉. 所以 i Row2 j Row1 i j Matrix f@ f* 放入 i Row2 (指令 i j i j matrix f@ Row2Row1Scalar*Substract>Row2)
這裡有經過測試, 方法一兩列都變成一後再相減消去會比較準確些, 所以這裡採用這樣的做法.
消去一列後變這樣
反覆操作後 Col: 0 下方都全部消去
同樣的程序, 透過外面的 do-loop 迴圈, 會逐步消去主軸下方的所有元素! 就這樣, 完成所有 LowerTrianglize 的程序!
再來就是消去上三角形, 主軸上方所有的的元素 UpperTrianglize
非常類似的做法, 不過這時候的主軸元素都已經是1了, 所以用方法二, 直接乘上非主軸的數值相等後相減消去.
: UpperTrianglize
MatSelected RowNum@ 1
do
0 i
do
i j MatSelected Matrix f@ f0<>
陣列元素已經為零的話, 什麼都不用做!
if
i j <>
if i j 2dup MatSelected Matrix f@
Row2Row1Scalar*Subtract>Row2
then
then
-1
+loop
loop
;
完成上三角化消去程序後, 就得到我們的答案了! 陣列的最後一欄就是我們要的解答!
主軸值的調整.
其實不調整主軸, 也可以得到解答的! 但是因為計算中有浮點運算的誤差, 當兩個相近的數字相減消去時這樣的誤差會擴散出來. 某些 ill condition 下, 最後會造成最終的誤差很大.
如果每一列的主軸元素, 是這一欄裡面數值最大的話, 這樣的誤差就可以被減小.
所以這裡要利用規則一: 陣列裡列(Row)跟列(Row)的交換, 不會影響解的位置或數值.
AdjustPivot 會利用這個規則, 在每一欄裡去尋找最大的元素, 找到的話就跟主軸列整列交換, 讓主軸元素永遠是同欄最大的那個, 來盡可能地減少浮點運算的誤差.
: AdjustPivot ( --)
MatSelected RowNum@ 0
do i i i MatSelected Matrix f@ fabs ( ColMax F: max) 留初值
( search for max pivot)
MatSelected RowNum@ i
do i j MatSelected Matrix f@ fabs ( ColMax F: max new)
fover fover f< 跟初值比較,比較大的話就留下來
if fswap drop i then
fdrop
loop
( swap with max pivot) 最後留下來的那個就是最大值囉!
fdrop i RowSwap 跟現有的主軸列交換!
loop
;
執行的狀況
完整程式列表
\
\ linear equations solver (Gaussian Elimination)
\
\ Frank Lin 2018.05.08
\
\ reset and change words searching orders, FORTH words first
FORTH only FORTH also
\ create new vocabulary for grouping and storing new words
vocabulary LinearSolver
\ set LinearSolver as current definition vocabulary
LinearSolver definitions
: DefMatrix
create ( RowNum ColNum --)
2dup , , * floats allot
;
: ColNum@ ( maddr -- ColNum)
@
;
: RowNum@ ( maddr -- RowNum)
cell+ @
;
: Matrix ( Row Col maddr -- addr')
>r r@ @ ( Row Col ColNum | maddr)
rot * ( Col ColNum*Row | maddr)
+ floats [ 2 cells ] literal +
r> +
;
: .matrix ( maddr --)
cr
dup RowNum@ 0
do dup ColNum@ 0
do dup j i rot Matrix f@ f. space
loop cr
loop drop
;
: EmptyMatrix ( maddr --)
dup @ ( maddr ColNum)
swap cell+ dup @ ( ColNum maddr' RowNum)
swap cell+ -rot * floats
erase
;
\ defer define for matrix
: DummyMatrix
0 value
;
\ Ex:
\ DummyMatrix MatSelected
\
\ 4 5 MatrixDef ExperimentData
\ ExperimentData to MatSelected
\
DummyMatrix MatSelected
: RowScalar* ( f.scalar row --)
MatSelected ColNum@ 0
do dup i MatSelected Matrix dup f@ fover f* f!
loop drop
;
: RowUnitizebyCol ( Row Col --)
2dup MatSelected Matrix f@ ( scalar) fdup f0<>
if 1e0 fswap f/ drop RowScalar*
else 2drop then
;
: Row2Row1Scalar*Subtract>Row2 ( Row2 Row1 f.scalar --)
MatSelected ColNum@ 0
do over i MatSelected Matrix dup f@ ( r2 r1 r2adr f.sc f.r2)
fover ( r2 r1 r2adr f.sc f.r2 f.sc)
over i MatSelected Matrix f@ f* ( r2 r1 r2adr f.sc f.r2 f.r1*sc)
f- f!
loop 2drop fdrop
;
: LowerTrianglize
MatSelected RowNum@ 0
do
MatSelected RowNum@ i
do
i j MatSelected Matrix f@ f0<>
if
i j RowUnitizeByCol
i j =
if j j MatSelected Matrix f@ f0= abort" Pivot is ZERO!!"
else i j 1e0 Row2Row1Scalar*Subtract>Row2 then
then
loop
loop
;
: UpperTrianglize
MatSelected RowNum@ 1
do
0 i
do
i j MatSelected Matrix f@ f0<>
if
i j <>
if i j 2dup MatSelected Matrix f@
Row2Row1Scalar*Subtract>Row2
then
then
-1
+loop
loop
;
: RowSwap ( row2 row1 --)
2dup = if 2drop exit then
MatSelected ColNum@ 0
do over i MatSelected Matrix dup f@ ( r2 r1 r2adr f.r2)
over i MatSelected Matrix dup f@ ( r2 r1 r2adr r1adr f.r2 f.r1)
fswap f! f!
loop 2drop
;
: AdjustPivot ( --)
MatSelected RowNum@ 0
do i i i MatSelected Matrix f@ fabs ( ColMax f.max)
( search for max pivot)
MatSelected RowNum@ i
do i j MatSelected Matrix f@ fabs ( ColMax f.max f.new)
fover fover f<
if fswap drop i then
fdrop
loop
( swap with max pivot)
fdrop i RowSwap
loop
;
: SolveLinearEQ ( m.addr --)
to MatSelected
AdjustPivot
LowerTrianglize
UpperTrianglize
;
\ add LinearSolver to searching-orders
LinearSolver also
\ create new vocabulary for matrix data
vocabulary MatrixData
\ set MatrixData as current definition vocabulary
MatrixData definitions
create test1
4 5 , ,
1.1e0 f, 7.2e0 f, 1.3e0 f, 1.4e0 f, 1.5e0 f,
2.1e0 f, 2.2e0 f, 7.3e0 f, 2.4e0 f, 4.5e0 f,
1.1e0 f, -3.2e0 f, 7.3e0 f, 3.5e0 f, -2.0e0 f,
4.1e0 f, 4.5e0 f, 5.1e0 f, -4.4e0 f, 7.5e0 f,
create test2 \ for ill condition
2 3 , ,
0.003e0 f, 591400e0 f, 591700e0 f,
5.291e0 f, -6.13e0 f, 46.78e0 f,
create test3
4 5 , ,
-2.0e0 f, 1.1e0 f, -2.0e0 f, -1.8e0 f, 1.0e0 f,
3.2e0 f, 2.1e0 f, 3.2e0 f, 2.2e0 f, 1.0e0 f,
3.4e0 f, 2.3e0 f, 4.1e0 f, 3.2e0 f, 6.0e0 f,
2.6e0 f, 1.1e0 f, -3.2e0 f, 2.4e0 f, -7.0e0 f,
xxx
留言列表
