2.5 多解矩陣方程的求解

盡管前面介紹的vpasolve()函數能夠一次性求出某些方程全部的解，但對一般矩陣方程，尤其是非線性矩陣方程而言是無能為力的，也沒有其他MATLAB程序能夠求解任意的非線性矩陣方程。

前面介紹了由給定初值求解函數的幾種方法，但這些方法很難一次性求解多解非線性方程所有的解，所以應該構建更簡單的函數完成這樣的任務，本節給出一種求解的思路，并依照該思路編寫出MATLAB通用程序，試圖得出方程感興趣區域內全部的解，并再擴展一步該方法，試圖得出方程的全部準解析解。

2.5.1 方程求解思路與一般求解函數

由前面給出的求解函數可見，如果選定了一個初值，則可以通過隨機數矩陣生成的方式產生一個初始搜索點，得出方程的解。更一般地，可以建立一個循環結構實現這樣的操作。由初始搜索點，調用fsolve()函數得出方程的一個解，如果這個解已經被記錄，則可以比較這個解和已記錄解的精度，如果新的解更精確，則用這個解取代已記錄的解，否則舍棄。如果這個解是新的解，則記錄該解。

如果這個循環結構設計成死循環，則有望得出方程全部的解。根據這樣的思路可以編寫出一個通用的求解函數，其實這個函數以前曾經發布了多個版本，這個版本中，特別增加了一些處理，例如，可以嘗試零矩陣是不是方程的孤立解；如果找到的根比以前存儲的更精確，則替換該根；如果找到的新解為復數，則檢驗其共軛復數是不是方程的根。基于這樣的考慮，編寫了下面的求解函數。這個函數有一個特點，運行的時間越長，得出的結果可能越精確。

該函數調用底層公用函數default_vals()，其內容在其他卷中有過描述，為方便起見，這里重新給出。

方程求解函數more_sols()的調用格式為

    more_sols(f,X0,a,?,tlim,opts)

式中，f為方程的函數句柄，可以由匿名函數與M函數描述原代數方程；X0為三維數組，用于描述解的初值，如果首次求解方程，建議將其設置為zeros(n,m,0)，即空白三維數組，n和m為解矩陣的維數；方程的解被自動存在MATLAB工作空間中的三維數組X中，如果想繼續搜索方程的解，則應該在X0的位置填寫X；a的默認值為1000，表示在[?500,500]區間內大范圍搜索方程的解；?的默認值為eps；tlim的默認值為30，表示30s沒有找到新的解就自動終止程序；還可以指定求解的控制選項opts，默認值為optimset。a還可以取為復數，表示需要求取方程的復數根。另外，a還可以給定為求解區間[a,b]。

例2-36 試重新求解例2-10中的一元超越方程。

解 從圖2-5給出的曲線可見，該方程在期望的區域內有6個交點，利用編寫的more_sols()函數可以求出這6個交點，并在圖2-11上直接標注出來。得出方程的解為x1=[1.4720,4.6349,7.7990,10.9601,14.1159,17.2666]。

例2-37 試求解例2-11給出的代數方程在?2π≤x1,x2≤2π區域內全部的根。

解 利用下面自然的方式可以直接求解非線性代數方程，先用匿名函數的形式描述聯立方程，這樣就可以調用more_sols()函數求取方程在感興趣區域內全部的數值解。這里將A可以選作13，比4π略大的值，得出的解個數多于在感興趣區域內的解。

可以提取出感興趣區域內全部的根，可以發現，總的根的個數為110個，可以將得到的解疊印到圖解法的圖形上，如圖2-12所示。

例2-38 試用數值方法重新求解例2-32中的廣義Riccati方程。

解 和以前一樣，先輸入幾個矩陣，然后由匿名函數的方式描述方程，再給出求解命令就可以直接求解方程，得到方程的8個實數根。

繼續求解方程則可以得出方程全部的復數根，如果將A設置成復數量。這樣總共可以得到方程的20個根，與例2-32得出的結果完全一致。

    >> more_sols(F,X,1000+1000i)

例2-39 如果例2-32中的廣義Riccati方程變換成DX+XA?BX2+C=0，試用不同的方法求解該方程。

解 如果用more_sols()函數直接求解，則可以得到19個實數根。

如果在復數范圍內搜索方程的根，總共可以找到57個根。現在嘗試vpasolve()函數的準解析解求解方法，經過5873.8s的等待，可以得到60個根。

兩種方法得出的根的數目差不多，不過仔細對比得出的根，如x11的值，可以發現，由后者得出的根有11個位于無窮遠處，其幅值大于1010，而前者得出的都是在合理范圍內的數值。兩種方法得出的x11對比在表2-1中給出。此外，經檢驗，x11=38的根是第21個根，將其提取并代入原方程，得出誤差矩陣的范數為18981334.363，不滿足原方程。x11=?1.5146的根是第15個根，也不滿足原始方程。由此可見，采用vpasolve()函數求解，即使處理多項式矩陣方程，得出的結果也是不可靠的。

例2-40 試求解非線性矩陣方程。

eAXsinBX?CX+D=0

其中A、B、C和D矩陣在例2-32中給出，試求出該方程的全部實根。

解 可以用下面的語句直接求解這里給出的復雜非線性矩陣方程，已經找到122個實根。用戶還可以自己嘗試，看看能不能得出更多的實根（根的存儲文件為data2_36.mat）。

2.5.2 偽多項式方程的求解

偽多項式方程是非線性矩陣方程的一個特例，因為未知矩陣實際上是一個標量。這里將給出偽多項式方程的定義，并給出求解方法。

定義2-10 偽多項式（pseudo-polynomial）方程的一般數學形式為

式中，αi為實數。

可見，偽多項式方程是常規多項式方程的擴展，求解方法可能遠難于普通多項式方程的求解方法。這里將探討不同方法的可行性。

例2-41 試求解偽多項式方程[5]x2.3+5x1.6+6x1.3?5x0.4+7=0。

解 一種容易想到的方法是引入新變量z=x0.1，這樣原方程可以映射成關于z的多項式方程如下：

z23+5z16+6z13?5z4+7=0

可以求出，該方程有23個根，再用x=z10就可以求出方程全部的根。這樣的思想可以由下面的MATLAB語句實現。

不過，把這樣得出的解代回到原來的方程可以發現，絕大部分的根都不滿足原有的偽多項式方程。原方程到底有多少個根呢？上面得到的x只有兩個根滿足原方程，即x=?0.1076±0.5562j，其余的21個根都是增根。由下面語句也可以得出同樣兩個根。

從數學角度看，這對真實的根位于第一Riemann葉上。其余的解位于其他Riemann葉上，都是方程的增根，但不滿足原方程。

例2-42 試求解無理階次偽多項式方程。

解 由于這里的階次是無理數，無法將其轉換成前面介紹的普通多項式方程，所以more_sols()函數就成了求解這類方程的唯一方法，可以由下面的語句直接求解該方程。可見，該無理階偽多項式方程只有兩個根，位置為s=?0.0812±0.2880j。

將得出的解代回原方程，則可見誤差為9.1551×10?16，該解在數值意義下足夠精確。從這個例子還可以看出，即使階次變成了無理數，并未給求解過程增加任何麻煩，求解過程和計算復雜度與前面的例子完全一致。

2.5.3 高精度求解函數

在這里給出的more_sols()函數中，核心求解工具是fsolve()函數，若將其替換為高精度的vpasolve()函數，則可以編寫出高精度的非線性函數求解程序。

該函數的調用格式為more_vpasols(f,X0,A,tlim)，其中還嵌入了為其設計的底層支持子函數sol2vec()，將得出的解轉換成行向量。輸入變量f可以為符號型的行向量來描述聯立方程，初始矩陣X0指定為zeros(0,n)，其中，n為未知數的個數。其他的輸入變元與前面介紹的more_sols()函數是一致的。返回的變量X(i,:)存儲找到的第i個解。值得指出的是，more_vdpsols()函數的運行速度比more_sols()函數慢得多，但精度也高得多。

例2-43 考慮例2-11中的聯立方程，試找出?2π≤x,y≤2π范圍內所有準解析解。

解 可以用下面的命令直接求解聯立方程：

要檢驗得出方程根的精度，則首先應該提取出感興趣區域內的根x0和y0，并對其進行排序，得出的根代入原方程后的誤差范數為7.79×10?32，比例2-37中得出的精度要高得多，所需的時間在半個小時左右，也遠高于more_sols()函數，用這樣的方法只找到區域內的105個根。得出的圖形類似于圖2-12，不過有幾個點缺失。

例2-44 試求解例2-40的高精度數值解。

解 可以試用下面語句求解方程，不過在用符號表達式描述方程時，并不能計算出方程的表達式f，所以不能調用后續的more_vpasols()函數，不能得出方程的高精度數值解，求解這類方程只能采用例2-40的普通數值解方法。