求解函数多项式型
⾮多项式型⼀维⾼维符号数值算法solve ⽀持，得到全部符号解若可符号解则得到根
⽀持⽀持⽀持当⽆符号解时 符号解⽅法：利⽤等式性质得到标准可解函数的⽅法基本即模拟⼈⼯运算vpasolve ⽀持，得到全部数值解（随机初值）得到⼀个实根
⽀持⽀持\times ⽀持未知fsolve 由初值得到⼀个实根由初值得到⼀个实根⽀持⽀持\times ⽀持
优化⽅法，即⽤优化⽅法求解函数距离零点最近，具体⽅法为信赖域⽅法。
包含预处理共轭梯度(PCG)、狗腿(dogleg)算法和Levenberg-Marquardt 算法等
fzero 由初值得到⼀个实根由初值得到⼀个实根⽀持\times \times ⽀持
⼀维解⾮线性⽅程⽅法，⼆分法、⼆次反插和割线法的混合运⽤
具体原理见数值求解⾮线性⽅程的和
roots ⽀持，得到全部数值解\times ⽀持\times \times ⽀持
特征值⽅法，即将多项式转化友矩阵(companion matrix)
然后使⽤求矩阵特征值的算法求得所有解（那是另外⼀个问题了）⾮线性⽅程（组）：MATLAB 内置函数solve,vpasolve,fsolve,fzer 。。。
MATLAB 函数 solve, vpasolve, fsolve, fzero, roots 功能和信息概览
也就是说，之前写了⼏篇关于⾮线性求解的，如⼆分法、⽜顿法（参见）、⼆次反插法（参见），只有⼀个库函数⽤了类似的⽅法。各函数⽤法详解
1. (符号/数值)解⽅程(组)函数 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
ans =
TRUE
TRUE
在solve中，还可以使⽤ assume 函数来规定符号变量的性质（如整数、⼤于零、区间等等）。
2. 多初值的数值解⽅程(组)函数 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
很轻松地得到了该⽅程的全部解⽽不⽤再去⼿动猜测了。
3. 数值解⽅程(组)函数 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)=0x_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 1
5 4.11698e-22 7.0962e-0
6 3.06e-11 1.61
输出内容中，iteration为迭代次数，func-count为函数的总调⽤次数，f(x)为函数值的⼀个性质（暂时还没搞清楚是啥，毕竟⼆维映射不可能只有⼀个值），Norm of step应当是迭代步长（相邻迭代点间隔）的范数（模长），first order optimality ⼀阶优化条件，最终迭代是否终⽌的判据就是⼀阶优化条件是否⾜够接近零。绘图可以看出，随着迭代的进⾏，⼀阶优化条件趋于零。
理论上，fsolve函数还允许指定求解的算法，⽐如使⽤单纯信赖域，或者使⽤狗腿信赖域，或者使⽤Levenberg-Marquardt算法。但总⽽⾔之，fsolve的算法均属优化算法，也因此在这⾥不⾜以讨论完全，留待优化部分的笔记说明。
4. 数值求⼀维函数零点函数 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;
matlab求导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
官⽅参考页：
除了求多项式根啥也⼲不了的⼀个函数，输⼊也和其他求根函数迥异。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
Processing math: 0%
发布评论