非線性方程(組)：MATLAB內置函數 solve, vpasolve, fsolve, fzero, roots [MATLAB]

MATLAB函數 solve, vpasolve, fsolve, fzero, roots 功能和信息概覽

求解函數	多項式型	非多項式型	一維	高維	符號	數值	算法
solve	支持，獲得所有符號解	若可符號解則獲得根	支持	支持	支持	當無符號解時	符號解方法：利用等式性質獲得標準可解函數的方法 基本即模擬人工運算
vpasolve	支持，獲得所有數值解	（隨機初值）獲得一個實根	支持	支持	$\times$	支持	未知
fsolve	由初值獲得一個實根	由初值獲得一個實根	支持	支持	$\times$	支持	優化方法，即用優化方法求解函數距離零點最近，具體方法爲信賴域方法。 包含預處理共軛梯度(PCG)、狗腿(dogleg)算法和Levenberg-Marquardt算法等
fzero	由初值獲得一個實根	由初值獲得一個實根	支持	$\times$	$\times$	支持	一維解非線性方程方法，二分法、二次反插和割線法的混合運用 具體原理見數值求解非線性方程的二分法、不動點迭代和牛頓法和插值方法
roots	支持，獲得所有數值解	$\times$	支持	$\times$	$\times$	支持	特徵值方法，即將多項式轉化友矩陣(companion matrix) 而後使用求矩陣特徵值的算法求得全部解（那是另一個問題了）

 

 　　也就是說，以前寫了幾篇關於非線性求解的，如二分法、牛頓法（參見二分法、不動點迭代、牛頓法）、二次反插法（參見插值法），只有一個庫函數用了相似的方法。

各函數用法詳解

1. (符號/數值)解方程(組)函數 solve

　　官方參考頁：Equations and systems solver - MATLAB solve

　　solve 是基本的用於符號解方程的內置函數，返回類型爲符號變量矩陣($m\times n$ sym)。當沒法符號求解時，拋出警告並輸出一個數值解。基本形式爲：

solve(eqn, var, Name, Val);　　% eqn爲符號表達式/符號變量/符號表達式的函數句柄, var爲未知量; Name爲附加要求，Val爲其值

　　能夠用solve解一維方程。對於多項式，solve能夠返回其全部值。

func1 = @(x)x^3 - 20*x^2 - 25*x + 500;    % 建立函數句柄.句柄中的變量不屬符號變量，不須要定義！
syms x exp1;    % 定義符號變量 x, exp1；
exp1 = x^3 - 20*x^2 - 25*x + 500;    % 符號表達式，包含符號變量. 符號變量必須先有上一行定義！

solve(exp1 == 0, x)    % 命令行輸入a，傳入一個包含符號表達式的等式，x爲所要求的變量
solve(exp1, x)    % 命令行輸入b，傳入一個符號表達式，函數默認求其零點
solve(func1(x), x)    % 命令行輸入c，傳入參數func1(x)等價於傳入了符號表達式，和輸入b徹底同樣
solve(func1(x) == 0, x)    % 命令行輸入d，這句話和a徹底同樣
solve(func1, x)    % 命令行輸入e，傳入參數func1，這是一個函數句柄，函數默認求其零點

ans =    % 命令行輸出，三個解，爲3*1的符號向量。對以上五種輸入輸出都徹底同樣
-5
5
20

　　對於不可符號求解的函數零點/方程解，solve拋出警告並返回一個數值解：

exp1 = atan(x) - x - 1;    % 不可符號求零點的表達式
solve(exp1 == 0, x)    % 命令行輸入
% 命令行輸出：
警告: Cannot solve symbolically. Returning a numeric approximation instead.
ans =
-2.132267725272885131625420696936

　　值得注意的是，雖然稱之爲「數值解」，而且輸出了一個位數很是多的小數，但查看數據類型發現ans的數據類型依然是符號變量。其實，若是是正常的double類型的變量，直接輸出是不可能給出這麼多位的。matlab裏面默認打印精度是4位小數，能夠用  format long  語句調整到15位小數，和機器精度基本持平。

　　solve也能夠求解方程組，此時輸入的表達式epn和變量var爲行向量（親測列向量不可用）：

exp1 = [x^2/9 + y^2/4 == 1, 2*x - 3*y == 0];    % 聯立橢圓方程和直線方程
solution = solve(exp1, [x, y]);    % 解方程組
% 命令行輸出
solution = 
    包含如下字段的struct:
    x: [2*1 sym]
    y: [2*1 sym]
% 這意味着x和y均有兩個解。函數輸出的solution是一個struct，訪問方法和C語言訪問struct成員同樣：
struct.x    % 命令行輸入
ans =    % 命令行輸出
 (3*2^(1/2))/2
-(3*2^(1/2))/2

　　能夠看出，當solve給出符號解的時候，它是不會化簡計算的。任何matlab的符號計算，包括四則運算、求導積分，都不具有化簡/數值計算的功能。

　　此外，solve函數還有一些選項可選，這使得符號求解名副其實，這纔是solve的強大之處。運用solve函數，Name設定爲'ReturnConditions'，其值設定爲true，表示要求solve函數輸出詳細信息。用這個方法咱們能夠得出一族解：

% 求解方程sin(x)=cos(x)，咱們知道這個方程有無窮多解
[solution, parameter, condition] = solve(sin(x)==cos(x), x, 'ReturnConditions', true);
% 命令行輸出
solution = pi/4 + pi*k    % 獲得一族解，以pi爲週期
parameter = k    % parameter輸出的是構成解的參數（符號變量）
condition =
in(k, 'integer');    % condition代表parameter的條件，此處k爲整數

　　而通常地，對於多變量的多項式（組），當多項式數量不足以肯定全部參數時，按照以上設定，solve函數能夠解出幾個變量關於其餘變量的函數：

exp1 = [x^2/9 + y^2/4 + z^2 == 1, 2*x - 3*y == 0];    % 橢球面和平面方程聯立，結果應爲曲線而非點
solution = solve(exp1, [x, y]，'ReturnConditions', true);    % 要求解出其中x和y的表達式
solution.x    % 命令行輸入：訪問solution結構體的x參數
ans =    % 命令行輸出：x關於z的表達式，是符號向量，能夠經過索引solution.x(1)和solution.x(2)分別訪問
 (3*2^(1/2)*(-(z - 1)*(z + 1))^(1/2))/2
-(3*2^(1/2)*(-(z - 1)*(z + 1))^(1/2))/2
% 結構體還有參數parameters和conditions。此處沒有新生成的參數，所以parameters和conditions沒有意義
solution.parameters    % 命令行輸入
ans = 
Empty sym: 1-by-0
solution.conditions    % 命令行輸入
ans = 
TRUE
TRUE

　　在solve中，還可使用  assume 函數來規定符號變量的性質（如整數、大於零、區間等等）。

 

2. 多初值的數值解方程(組)函數 vpasolve

　　官方參考頁：Solve equations numerically - MATLAB vpasolve

　　vpasolve是一款數值解方程(組)的函數。輸入一些個參數，返回符號型數值標量/向量（即以數值的形式顯示但實際上仍是符號變量）。它的基本使用方式是：

vpasolve(eqns, vars, init_guess, 'Random', randomvalue);    % 方程(組)eqns，變量vars，初值點init_guess(可缺省，在random模式下可寫區間)，'Random'設置randomvalue（可缺省）

　　它的輸入、功能和輸出都和solve相仿。方程組的輸入一樣爲行向量，變量組的輸入也同樣。

　　當輸入一個能夠定解的多項式方程（組）時，vpasolve將會直接給出方程的數值解；若多項式方程數量不足以肯定全部的解，那麼vpasolve將會給出以剩餘變量表示的所求變量的函數，只是表達式的一部分（係數等）可能會以數值的形式呈現。注意，有理分式方程將會多項式化之後同樣處理。對於這些方程，init_guess的值寫了也沒用。

exp1 = (x-1)*(x-2)/(x-3);    % 分式方程
solution = vpasolve(exp1, x)    % 命令行輸入
solution =    % 命令行輸出，獲得了所有解
1.0
2.0

　　對於多元方程組，vpasolve的輸出也是struct結構體，訪問方法也和solve輸出的struct同樣。不一樣的是，vpasolve沒法求出含參的解，即沒法設定'ReturnConditions'選項。和solve相似，除了多項式方程和有理分式方程之外的任何方程，vpasolve都不會給出所有解。這樣一看，彷佛vpasolve只不過就是把solve的結果所有轉化爲數值形式，甚至許多solve的功能都不能知足。這樣的想法固然不對，vpasolve也有其自身的優點，這來源於：

　　A）能夠設置初值進行數值求解。對於不可符號求解的方程，solve由於沒有設定初值，可能沒法搜索到合適的解。vpasolve則能夠設置初值，從而能夠進行後續解的搜索；B）能夠隨機取初值。咱們都知道求解方程和最優化問題的初值選取很是玄學，而一樣的初值最多隻能有一個解。而結合循環等控制語句，利用vpasolve的隨機初值功能可讓這一過程變得比較簡單。好比能夠寫做初等函數的半整數階Bessel(貝塞爾)函數，其零點有無窮多，但沒法經過符號方法求解，在solve中會遇到很大的問題，可是用vpasolve設置合適的初值能夠獲得許多組臨近初值的解。好比：

syms x;
exp1 = sin(x)/x;
exp1 = diff(exp1, x);    % 對原函數求導，獲得的函數零點和3/2階貝塞爾函數的零點相同
% 命令行輸入
solution = solve(exp1, x)
% 命令行輸出，每次獲得的結果均同樣，爲一個很遙遠的解
警告: Cannot solve symbolically. Returning a numeric approximation instead.
solution = 
-227.76107684764829218924973598808

solution = vpasolve(exp1, x, 1)    % 命令行輸入
solution =     % 命令行輸出大概就是0
0.00000000000000000000000014187233415524662526214503949551

solution = vpasolve(exp1, x, 3)    % 命令行輸入
solution =     % 命令行輸出
4.4934094579090641753078809272803　　

 　　另外，一些有限個解的方程，好比 $atanx=x/2$ ，咱們已經知道它有解0，根據畫圖還能夠肯定在x>0和x<0範圍內各有一個解。根據atanx的性質，咱們知道全部的解確定在區間[-5,5]之中。若是使用solve，每次均有警告而且輸出同樣，沒法得到三個不一樣的解；即便是以後講的fsolve也須要每次給定初始估計(init guess）。對於vpasolve，當肯定範圍了之後能夠簡單地使用循環的控制語句，只須要規定隨機撒點的區間爲[-5,5]：

syms x;
exp1 = atan(x) - x/2;
for i = 1:5
    solution = vpasolve(exp1 == 0, x, [-5, 5], 'Random', true);
    disp(solution);
end

　　輸出結果：

-1.3628993400364241574879337535051e-53    % 也差很少即0
2.3311223704144226136678359559171    % 這就是要求的正根
-2.3311223704144226136678359559171    % 這就是要求的負根，和正根關於原點對稱
2.3311223704144226136678359559171
0

　　很輕鬆地獲得了該方程的所有解而不用再去手動猜想了。

 

3. 數值解方程(組)函數 fsolve

　　官方參考頁：Solve system of nonlinear equations -  MATLAB fsolve

 　　fsolve多是目前matlab的內置庫函數中最經常使用的求（非線性）方程（組）解的函數，也是最爲人熟知的。它用於數值求解方程（組），具備較廣的適用範圍（適用於高維和非線性、非多項式情形），甚至能夠求矩陣方程的解（即甚至能夠求解未知量爲矩陣的情景）。fsolve函數的基本形式是：

[x, fval, exitflag, output, jacobian] = fsolve(func, init_guess, options); 
% 輸入函數句柄func，初值（向量）init_guess，求解選項options（可缺省）;
% 輸出解x，x處值fval（也就是殘差，可缺省），計算終止信息exitflag、輸出信息output、x處雅可比矩陣jacobian（都可缺省）

　　好比參考頁面給出的示例非線性方程組：$$e^{-e^{x_1+x_2}}-x_2(1+x_1^2)=0$$ $$x_1cos(x_2)+x_2sin(x_1)=\frac{1}{2}$$　　這是一個迷通常的方程組，嵌套的天然指數讓人十分混亂，咱們也並不指望獲得這個方程的符號解或者解析解。咱們將該方程組轉化爲matlab函數句柄：

x = sym('x', [1,2]);
% 生成符號變量向量
f = [exp(-exp(-x(1)+x(2))) - x(2)*(1 + x(1)^2), x(1)*cos(x(2)) + x(2)*sin(x(1)) - 0.5];
% 生成函數句柄func，該句柄的輸入參數爲一貫量
func = matlabFunction(f, 'Vars',{[x(1), x(2)]});

　　而後調用fsolve對於函數func進行求解，輸出一個求解消息和解solution：

% 命令行輸入
solution = fsolve(func, [0,0])
% 命令行輸出
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.
solution =
    0.3931    0.3366

　　須要注意的是，fsolve輸入的函數句柄func只接受一個變量！fsolve可用於高維的情形，如例子中的二維，是經過將函數句柄的輸入轉化爲向量實現的，即func接受一個向量形式的變量。對於建立一個輸入參數爲向量的函數句柄，簡單地採用@方法彷佛是行不通的。以上採用的方法是利用函數matlabFunction，定義變量('Var')爲向量[x(1),x(2)]，從而將符號變量f轉化爲函數句柄func。另外一種可能更加普適（但更加麻煩）的方法參見官方參考頁的示例或者matlab中函數fsolve的help文檔，經過定義一個函數文件來實現這一操做（函數function文件和函數句柄是等價的）。

　　函數fsolve提供了一些能夠做爲輸出設置的options選項。options的設置以下：

options = optimoptions('fsolve', 'Display', opt1, 'PlotFcn', opt2);
% opt1能夠填 'none'/'off'（無額外顯示）/'iter'（迭代信息表格）
% opt2能夠填函數 @optimplotfirstorderopt 繪製一階極值條件隨迭代的變化

　　如今，嘗試使用'iter'和'@optimplotfirstorderopt選項：

options = optimoptions('fsolve', 'Display', 'iter', 'PlotFcn', @optimplotfirstorderopt);
solution = fsolve(func, [0,0], options);
% 除了上述輸出，還有了另外的輸出：
                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          3        0.385335                         0.503               1
     1          6       0.0336737       0.642449          0.206               1
     2          9     0.000110235       0.122659         0.0162            1.61
     3         12     8.13142e-11     0.00681475       1.13e-05            1.61
     4         15     4.11698e-22     7.0962e-06       3.06e-11            1.61

　　

　　輸出內容中，iteration爲迭代次數，func-count爲函數的總調用次數，f(x)爲函數值的一個性質（暫時還沒搞清楚是啥，畢竟二維映射不可能只有一個值），Norm of step應當是迭代步長（相鄰迭代點間隔）的範數（模長），first order optimality 一階優化條件，最終迭代是否終止的判據就是一階優化條件是否足夠接近零。繪圖能夠看出，隨着迭代的進行，一階優化條件趨於零。

　　理論上，fsolve函數還容許指定求解的算法，好比使用單純信賴域，或者使用狗腿信賴域，或者使用Levenberg-Marquardt算法。但總而言之，fsolve的算法均屬優化算法，也所以在這裏不足以討論徹底，留待優化部分的筆記說明。

 

4. 數值求一維函數零點函數 fzero

　　官方參考頁：非線性函數的根 - MATLAB fzero

　　fzero用於求函數零點。這個函數致力於求解一維函數的零點。其基本形式：

x = fzero(func, init_guess, options)    % func爲函數句柄，init_guess爲初值，options能夠包括其餘要求（可缺省）

　　fzero在應用上最使人高興的是其豐富的輸出內容，有利於觀察迭代的結果，這用到options控制。options的控制方法爲：

options = optimset(A, B);    % A爲一個選項名，B爲該選項值

　　而後將options變量帶入函數便可。具體能夠參見參考頁，在此舉兩個例子，好比但願輸出迭代的每一步：

options = optimset('Display', 'iter');　　% 設定'Display'選項爲'iter'模式
func = @(x)x^3 + x + 5;
zeropt = fzero(func, 10, options);
disp(zeropt);

　　則有輸出（節選）：

Func-count    x          f(x)             Procedure
   26        -4.05097      -65.5287        initial
   27        -3.40338      -37.8248        interpolation
   28          -2.541      -13.9473        interpolation
   29          -2.541      -13.9473        bisection
   30        -1.65947      -1.22938        interpolation
   31        -1.57533     -0.484774        interpolation
   32        -1.52086    -0.0386585        interpolation
   33        -1.51616   -0.00138072        interpolation
   34        -1.51598  -4.15907e-06        interpolation
   35        -1.51598  -4.49535e-10        interpolation
   36        -1.51598  -8.88178e-16        interpolation
   37        -1.51598  -8.88178e-16        interpolation

　　從中，咱們能夠看到每一步的x變化，f(x)的取值，甚至每一次迭代執行的操做：是二分法(bisection)或者插值類方法(interpolation)。咱們還能夠將迭代步驟可視化：

options = optimset('PlotFcns', @optimplotfval);　　% 每次輸出函數值
func = @(x)x^3 + x + 5;
zeropt = fzero(func, 0, options);
disp(zeropt);

　　輸出圖片：

 

5. 數值求多項式零點函數 roots

　　官方參考頁：多項式根 - MATLAB roots

　　除了求多項式根啥也幹不了的一個函數，輸入也和其餘求根函數迥異。roots的標準形式以下，輸入一個行向量，返回一個double型的列向量：

r = roots(p);    % 其中，p爲一個行向量，裏面依次爲多項式降次排序時各次項係數（若無該次則寫0）

 

　　roots也不是一無可取。相比於fzero和fsolve這樣的函數，roots能夠給出多項式的全部解，包括實數解和複數解：

p = roots([1, 0, 0, -1])    % 命令行輸入
p =     % 命令行輸出三個解，其中一對共軛復根，一個實根
  -0.5000 + 0.8660i
  -0.5000 - 0.8660i
   1.0000 + 0.0000i