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