MATLAB函数 solve, vpasolve, fsolve, fzero, roots 功能和信息概览html
求解函数 | 多项式型 | 非多项式型 | 一维 | 高维 | 符号 | 数值 | 算法 |
solve | 支持,获得所有符号解 | 若可符号解则获得根 | 支持 | 支持 | 支持 | 当无符号解时 | 符号解方法:利用等式性质获得标准可解函数的方法算法 基本即模拟人工运算app |
vpasolve | 支持,获得所有数值解 | (随机初值)获得一个实根 | 支持 | 支持 | $\times$ | 支持 | 未知 |
fsolve | 由初值获得一个实根 | 由初值获得一个实根 | 支持 | 支持 | $\times$ | 支持 | 优化方法,即用优化方法求解函数距离零点最近,具体方法为信赖域方法。dom 包含预处理共轭梯度(PCG)、狗腿(dogleg)算法和Levenberg-Marquardt算法等函数 |
fzero | 由初值获得一个实根 | 由初值获得一个实根 | 支持 | $\times$ | $\times$ | 支持 | 一维解非线性方程方法,二分法、二次反插和割线法的混合运用优化 具体原理见数值求解非线性方程的二分法、不动点迭代和牛顿法和插值方法 spa |
roots | 支持,获得所有数值解 | $\times$ | 支持 | $\times$ | $\times$ | 支持 | 特征值方法,即将多项式转化友矩阵(companion matrix)命令行 而后使用求矩阵特征值的算法求得全部解(那是另一个问题了)code |
也就是说,以前写了几篇关于非线性求解的,如二分法、牛顿法(参见二分法、不动点迭代、牛顿法)、二次反插法(参见插值法),只有一个库函数用了相似的方法。orm
各函数用法详解
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