現代營建 454---Sylvester矩陣方程式求解：(使用Windows Presentation Foundation的執行環境)

Sylvester矩陣方程式求解：

(使用Windows Presentation Foundation的執行環境)

□葉明養

一、前言

        有很多的學問和研究，從表面看起來似乎平凡無奇，但是如果真正深入的去研究和學習，終其人生都無法窺究其全貌，學習研究有所謂deep-first和breadth-first等的學習方式，兩者何者較優的學習，其實是一體的兩面。對個人而言，人的一生有限，學習研究並精通某一學問，所花費的時間和學問知識的了解，就像是無底洞一樣，永遠達不到底，如何有效的利用學習，也是人生的另一項難題。
    單單就C#程式語言一項而言，要百分百完全掌控、能完全應用自如都不是一件容易的事情，我個人認為甚至C#的建構者(Creater)，Anders Hejlsberg先生也都無法百分百完全掌控C#程式語言，因為程式語言與作業系統相互配合，並需廣泛提供各式各樣的程式庫供使用者使用，語言編譯器(如現今微軟Roslyn分析器)需不斷的提升改進, 還有其它種種事情，都需要一群專家共同投入研發，更需要大財團的財力支援等等，並不是個人單獨能支撐開發完成的，無論是程式語言的開發者或是使用者，終其一生不可能完全了解且完全掌控應用自如。人的一生有限，但程式語言的生命週期循環是無限的，若不能被大眾所接受，或不能適應就被淘汰了。
      就我個人而言，民國六十二年，我選過虞兆中教授的高等結構學，林聰悟老師以Fortran程式語言教授矩陣結構，這是我學習結構矩陣和Fortran程式語言撰寫的啟蒙教育，但畢業就職後，從此就中斷也沒有機會再使用Fortran程式撰寫，那時候中大型計算機是學校與大型顧問公司才擁有，直到民國七十八年我才在台大電機系的計算機訓練中心學習，以Pascal程式語言撰寫及學習資料結構開始，直到今天才有一點點C#程式語言撰寫的基礎，前後都超過二十八年了。我們常常聽到結構勁度矩陣這個名詞，但我不知有幾個人真正使用程式語言撰寫建築結構物的勁度矩陣？大部分的人都使用Sap2000等相關的Package作結構分析, 而結構軟體僅能提供計算結果，不太可能提供中間過程數據的結構勁度矩陣，所以學習線性代數與矩陣計算一直是我希望的目標，雖然時至今天從未曾再接觸過Fortran程式語言，但使用C#程式語言撰寫結構計算和鋼結構桿件設計，我覺得更方便更容易，物件導向(Object-Oriented)比結構性(Structured)的程式語言更容易撰寫，而結構性比非結構性(Non-Structured)的程式語言更容易撰寫。
       茲將學習C#程式語言和學習線性代數與矩陣計算的小小心得報告，撰寫“Sylvester矩陣方程式的求解”乙文，並結合微軟Windows Presentation Foundation視窗環境的展現，投稿現代營建雜誌社。軟體的開發並非是電機資訊界的專利，希望土木結構工程界和學術界先進，大家互相共同努力，能為土木結構工程，再開創工程軟體的領域知識(Domain Knowledge)。

二、熟悉程式語言與熟練開發工具

       Matlab，Maple，和Mathematica是大家比較熟悉的商業數學軟體，這三種Package有其本身的程式語言，在網路上也常常看到自由軟體SageMath，其它如Netlib，Numerical Algorithms Group(NAG)，和Internal Mathematics and Statics Library(IMSL)等商業軟體，使用Fortran，C，Java，和C#程式語言撰寫數學程式庫。而我使用C#撰寫矩陣計算程式庫，並使用Matlab作必要的驗證。
       撰寫App須要對程式語言有相當程度的了解，以降低解決問題的困難度，否則一面要處理問題，又有語言上的問題，就像新手駕駛一樣，一面路況不熟，一面駕駛技術不好，會有更嚴重的挫折感，所以對語言的熟悉度是撰寫App的首要條件。
       如果有人想要自己使用程式語言，撰寫程式碼解決問題，到底使用何種語言較合適，當然是使用自己最熟悉的語言，我自己過去曾經使用Pascal，C，C++，和C#，就我個人而言，Pascal太舊而C太低階都不太適合撰寫矩陣的App，若使用微軟的Managed C++，倒不如使用C#更合適，C#是最適合撰寫線性代數與矩陣計算的程式語言，或許大家會懷疑，最主要的原因是矩陣是C#程式語言內建(Built-In)的型態，也是陣列(Array)的一種，以double[,]的方式來表示，矩陣的加、減、乘、求逆矩陣、各種矩陣的分解等等，就像一般浮點數的運算一樣，而一般Column Vector或是Row Vector也是當作矩陣來處理。
       舉一個線性代數的書中，最常看到簡單的例子A X = B，左式中A和B是已知的矩陣，則未知的矩陣X = Ai B，Ai是A的逆矩陣，Ai矩陣乘B矩陣等於X矩陣，如果 A x = b ，b是已知Column Vector，x是未知Column Vector，則x = Ai b，與上述的求解完全相同。如果懂得C#程式語言，撰寫矩陣計算App，變得非常容易了，即使不懂C#程式語言，也很容易求得未知矩陣X或是Column Vector x。就像Matlab矩陣計算軟體一樣，直接代入已知A和B矩陣，很容易求取矩陣X或Column Vector x了。
        C#程式的開發工具就是微軟Visual Studio(VS)，可直接從網路下載直接安裝使用，我並未使用其他語言的開發工具，所以無法比較優缺點。使用VS撰寫App時，若語法有錯誤，無需編譯即可顯示在螢幕上，非常方便，但如果是邏輯錯誤，VS雖然提供Debug的功能，但我並沒有使用此項功能，這種Debug太低階了，而我是把可能有錯誤點的地方印出，直接查看前後程式碼的關係，找出錯誤的地方，這種Debug的方式，到目前為止運行的很好，所以我並沒有使用VS提供的Debug的功能。

三、以CSharp撰寫和求解矩陣方程式經驗談

        Sylvester和Lyapunoc矩陣方程式是用於控制系統工程，我完全不懂控制系統，當初學習線性代數時，在有關矩陣計算的文章中看到⊗這個符號，以及Sylvester矩陣方程的求解中，都使用這個符號來表示，感覺很新奇。雖然控制系工程並不是我要學習的東西，但仍然嘗試著使用C#程式語言撰寫求解，並在參考文獻[4]中得到驗證。
        先介紹二個定義，第一是⊗表示Kronecker Product，若設定已知A是mxn的實數矩陣，B是pxq的實數矩陣，則A⊗B的定義是

       第二是Vec表示Vectorization，C 是mxn的實數矩陣，則
       C = [ C1 C2 . . . Cn ] 而 C1，C2， . . . Cn 都是m列的Column Vector
       Vec( C ) = [ C1 C2   . . . Cn ]_t 是mn x 1 的Column Vector。
       Vect(C) = [c₁₁, c₂₁, c₃₁, . . . cm₁，c₁₂, c₂₂, c₃₂，. . . ，cm₂，. . . ，c_1n， . . . ，c_mn]_t        _{是mn x 1的Column
Vector。

       有了A⊗B
和 Vec(C) 的觀念，就容易解釋Sylvester矩陣方程式了。

       Sylvester矩陣方程式是

       A X + X B = C                        . .
. . . . . . . . . . . ...............
( 1 )

       A X + X B

就是
A * X + X * B，其中A矩陣是mxm

，B矩陣是
nxn的正方形矩陣，C矩陣是
mxn的實數矩陣，X是待求取的未知矩陣。

       X未知的矩陣的求解方程式是

       ( I_b⊗A + B^t⊗I_a ) Vec( X ) = Vec ( C )
. . . . . . . . . . . ( 2 )

對程式開發者(Developer)而言，所謂演算法(Algorithms)，就是依據公式(2)分解成數個步驟，求得未知矩陣X，其求解步驟如下：

       步驟[1]：計算U1
= I_b}_{⊗A 。

       步驟[2]：計算U2
= B^t⊗I_a

(B^t

是
B 的
Transpose)。

       步驟[3]：計算U3
= U1 + U2，再設定U3A
= [U3]i ( Inverse of [U3] )。

       步驟[4]：計算U4
= Vec (C ) 。

       步驟[5]：計算U5
= U3A * U4 = Vec ( X )

。}
       步驟[6]：由Vec ( X )再求得X未知矩陣。
       公式(2)的求解方程式是數學家推導的結果，我不會而且也沒必要去推導，我是土木結構工程師，能看懂公式的意思，能使用CSharp程式語言去撰寫程式碼求解即可，公式求證與推導並不是我的專業，那就讓數學家們去做好了，就像我僅會駕馭車輛，但引擎蓋下(Under the Hood)的複雜機械與電機原理，實在不必完全去理解，就委由引擎專家去處理設計吧。
       我當初也懷疑公式(1)與公式(2)是否相同，但經撰寫C#程式碼，並從書本及參考文獻[4]的實例，測試驗證結果都完全符合，從公式(2)求取X矩陣較難，但從公式(1)驗證則非常容易，只要將X矩陣代入公式(1)的左邊，是否等於C矩陣，即可判斷對或錯了。
       其他相關的矩陣方程式公式 ( 3 )，相同的表示式公式 ( 4 )。
        A X B = C                             ................… . . . . . . . . . . . .( 3 )
       ( B^t⊗A ) Vec ( X ) = Vec ( C )             . . . . . . . . . . . . . . ( 4 )
       換言之，公式 ( 3 )之X矩陣，可由公式 ( 4 )求得。
       若A和B是正方形矩陣，且DET( A )和DET( B )不為0，則 X = A_i C B_i ，左式之i下標表示該矩陣的逆矩陣(Inverse Matrix)。
       若A是mxn矩陣，且m>n，B是p x q且q > p，則由廣義的逆矩陣可得。
       X = ( A^tA )_i ( A^tC B^t ) ( B B^t)_i ，但左式我並沒有使用C#程式去測試，也不知是否正確，希望將來有機會再撰寫程式碼求證，公式中之t上標表示該矩陣的轉置矩陣(Transpose Matrix)。
       Lyapunov矩陣方程式：
       A X + X A^t = C . . . . . . . . . . . . . . . . . . . . . . . . . . . .. ( 5 )
   公式( 5 )與公式( 1 )相同，僅設定 B = A^t 即可。

四、使用微軟視窗執行環境

使用C#程式語言求解Sylvester矩陣方程式，都是撰寫成Class的形式，並加入類別程式庫中，再編譯成實體Assembly DLL的檔案形式，而邏輯的形式就稱名稱空間(namespace)，直接加入應用程式中(using namespace)，一般較單純的執行環境是螢幕應用程式 (Console App)，這也是最普遍的執行環境，除了使用Console App的執行模式外，另外再撰寫視窗(Windows Presentation Foundation)的執行模式，無論是螢幕和視窗的執行模式，兩者都引入相同程式類別庫，顯然視窗的執行環境，能展現較豐富的使用者介面，但須要撰寫更多的程式碼。
在視窗程式中，視窗的左側是Sylvester Title和漸層色彩背景，展現較精緻的外表，右側是一個圖像，並在兩者之間加入一個GridSplitter，可以左右移動圖像，底部再加入兩個按鈕，分別顯示出三個已知矩陣和輸出計算結果矩陣，事實上視窗的展現可以有不同的變化，作者僅是Windows視窗的的初學者，仍待學習。最後在本文後段列出相關的程式碼提供參考。

五、結論

       許許多多自然界問題的求解，都把問題寫成矩陣的形式，然後求解矩陣就可得到解答，譬如多元方程式、常微分方程式、偏微分方程式、運動方程式等等求解，學習線性代數和矩陣計算是必備的課程，了解矩陣的解法，一般都到此為止就結束了。事實上真正的去求解矩陣又是另外一回事，求解矩陣不可能再使用手算，或許會使用Matlab數學軟體求解，上述Sylvester矩陣方程式是作者使用C#程式語言求解的例子。也藉此機會拋磚引玉，自己撰寫程式碼，求取適合自己的需求，而Matlab等等的Package並非能完全適合個人的需求。

輸出、視窗展現、螢幕App、和C#程式碼
(1)螢幕程式輸出 :

(2)螢幕程式碼：

using System;

using ClassLibrary1;

using MatData;

namespace T01

{

class Program

{

static void Main(string[] args)

{

// 已知矩陣[A]，[B]，[C]。

/* 以下是第一例的已知矩陣。

double[,] A = { {1, -1, 1}, {1, 1, -1}, {1, 1, 1} };

double[,] B = { { 8, 1, 6 }, { 3, 5, 7 }, { 4, 9, 2 } };

double[,] C = { { 1, 0, 0}, { 0, 1, 0}, { 0, 0, 1} };

// 以下是第二例的已知矩陣。

double[,] A = { {1, 0, 2, 3}, {4, 1, 0, 2}, {0, 5, 5, 6}, {1, 7, 9, 0} };

double[,] B = { {0, -1}, {1, 0} };

double[,] C = { {1, 0}, {2, 0}, {0, 3}, {1, 1} };

// Sylvester矩陣方程式計算。

// A X + X B = C ; A，B，C是已知矩陣，X是未知矩陣。

Sylvester syl = new Sylvester(A, B, C);

double[,] X = syl.Matrix;

// 將結果印在螢幕上。

Console.Write("\n** 矩陣[X] : **\n\n");

Console.Write(new PR(X));

}

** 矩陣[X] : **

0.12231 -0.07247 0.01311

-0.08057 -0.01615 0.15870

-0.01643 0.17835 -0.10723

請按任意鍵繼續 . . .

** 矩陣[X] : **

0.47316 -0.36643

-0.40057 0.35311

0.33053 -0.11422

0.07740 0.35601

請按任意鍵繼續 . . .

(3)視窗程式外觀：

(4)視窗程式輸出 :

已知A，B，和C三個矩陣 :

(Sylvester矩陣輸出及視窗加入GridSplitter可以調整左右圖示的大小)

(5)視窗程式碼 :

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.Threading.Tasks;

using System.Windows;

using System.Windows.Controls;

using System.Windows.Data;

using System.Windows.Documents;

using System.Windows.Input;

using System.Windows.Media;

using System.Windows.Media.Imaging;

using System.Windows.Navigation;

using System.Windows.Shapes;

using ClassLibrary1;

using MatData;

namespace W01

{

/// <summary>

/// Interaction logic for MainWindow.xaml

/// </summary>

public partial class MainWindow : Window

{

private double[,] A, B, C;

public MainWindow()

{

InitializeComponent();

// 已知[A],[B],[C]三個矩陣。

A = new double[,]

{ {1, 0, 2, 3}, {4, 1, 0, 2}, {0, 5, 5, 6}, {1, 7, 9, 0} };

B = new double[,]{ {0, -1}, {1, 0} };

C = new double[,]{ {1, 0}, {2, 0}, {0, 3}, {1, 1}};

// 視窗內容設定。

Title = "矩陣方程式求解";

Width = 750;

Height = 500;

// 建構主Grid(包含3個Column)，並加入視窗中。

Grid grid = new Grid();

for(int i = 0; i != 3; i++)

{

grid.ColumnDefinitions.Add(new ColumnDefinition());

}

grid.ColumnDefinitions[0].Width = new GridLength(500);

grid.ColumnDefinitions[1].Width =

new GridLength(0, GridUnitType.Auto);

Content = grid;

// 建構第二個Grid(包含2個Row)，並加入主Grid[0, 0]中。

Grid grid2 = new Grid();

for(int i = 0; i != 2; i++)

{

grid2.RowDefinitions.Add(new RowDefinition());

}

grid2.RowDefinitions[0].Height = new GridLength(400);

grid2.RowDefinitions[1].Height = new GridLength(50);

grid.Children.Add(grid2);

Grid.SetRow(grid2, 0);

Grid.SetColumn(grid2, 0);

// 建構TextBlock，並加入Grid[0, 0]=>Grid[0, 0]中。

TextBlock txtBlock = new TextBlock();

txtBlock.Text = "\nSylvester矩陣方程式求解:";

txtBlock.FontSize = 60;

txtBlock.FontWeight = FontWeights.Bold;

txtBlock.TextWrapping = TextWrapping.Wrap;

LinearGradientBrush brush = new LinearGradientBrush();

brush.GradientStops.Add(new GradientStop(Colors.Violet, 0));

brush.GradientStops.Add(new GradientStop(Colors.Yellow, 0.5));

brush.GradientStops.Add(new GradientStop(Colors.LightGreen, 0.8));

brush.GradientStops.Add(new GradientStop(Colors.Red, 1));

txtBlock.Background = brush;

txtBlock.HorizontalAlignment = HorizontalAlignment.Stretch;

txtBlock.VerticalAlignment = VerticalAlignment.Stretch;

grid2.Children.Add(txtBlock);

Grid.SetRow(txtBlock, 0);

Grid.SetColumn(txtBlock, 0);

// 建構StackPanel，並加入Grid[0, 0]=>Grid[0, 1]中。

StackPanel stack = new StackPanel();

stack.Orientation = Orientation.Horizontal;

stack.Margin = new Thickness(0);

brush = new LinearGradientBrush(

Colors.Red, Colors.Green, new Point(0, 0), new Point(1, 0));

stack.Background = brush;

//new SolidColorBrush(Colors.DarkOliveGreen);

grid2.Children.Add(stack);

Grid.SetRow(stack, 1);

Grid.SetColumn(stack, 0);

Button btn1 = new Button();

btn1.Content = "已知矩陣";

btn1.Margin = new Thickness(5);

btn1.Click += Btn1_Click;

stack.Children.Add(btn1);

Button btn2 = new Button();

btn2.Content = "計算結果";

btn2.Margin = new Thickness(5);

btn2.Click += Btn2_Click;

stack.Children.Add(btn2);

// 建構GridSplitter，並加入Grid[0, 1]中。

GridSplitter split = new GridSplitter();

split.ShowsPreview = true;

split.HorizontalAlignment = HorizontalAlignment.Center;

split.VerticalAlignment = VerticalAlignment.Stretch;

split.Width = 6;

grid.Children.Add(split);

Grid.SetRow(split, 0);

Grid.SetColumn(split, 1);

// 建構Button，並加入Grid[0, 2]中。

Uri uri = new Uri("file://C:\\2016\\Misc_05\\IMG_01.jpg");

BitmapImage bitmap = new BitmapImage(uri);

Image img = new Image();

img.Source = bitmap;

img.Stretch = Stretch.Fill;

Button btn = new Button();

btn.Content = img;

grid.Children.Add(btn);

Grid.SetRow(btn, 0);

Grid.SetColumn(btn, 2);

}

// 輸出計算結果。

private void Btn2_Click(object sender, RoutedEventArgs e)

{

string str = "Sylvester矩陣:\n";

Sylvester syl = new Sylvester(A, B, C);

double[,] X = syl.Matrix;

PR pr = new PR(X);

str += pr.ToString();

MessageBox.Show(str);

}

// 輸出三個已知矩陣。

private void Btn1_Click(object sender, RoutedEventArgs e)

{

string str = "[A]矩陣:\n";

PR pr = new PR(A);

str += pr.ToString() + "\n[B]矩陣:\n";

pr = new PR(B);

str += pr.ToString() + "\n[C]矩陣:\n";

pr = new PR(C);

str += pr.ToString();

MessageBox.Show(str);

}