非线性方程(组):MATLAB内置函数 solve, vpasolve, fsolve, fzero, roots [MATLAB]

摘要:
MATLAB函数求解,vpsolve,fsolve,fzero,根函数和信息概述求解函数多项式型非多项式型一维高维符号数值算法求解支持,获得所有符号解如果解可以签名,当没有符号解时获得根支持符号解方法:利用方程的性质获得标准可解函数的方法基本上是模拟手动操作vpsolve支持,获取所有数值解以获得实根支持$imes$support未知fsolve从初始值获取实根从初始值获得实根support$imes$support优化方法,即使用优化方法来求解最接近零点的函数,具体方法是信赖域方法。

MATLAB函数 solve, vpasolve, fsolve, fzero, roots 功能和信息概览

求解函数多项式型非多项式型一维高维符号数值算法
solve支持,得到全部符号解若可符号解则得到根支持支持支持当无符号解时

 符号解方法:利用等式性质得到标准可解函数的方法

基本即模拟人工运算

vpasolve支持,得到全部数值解(随机初值)得到一个实根支持支持$ imes$支持未知
fsolve由初值得到一个实根由初值得到一个实根支持支持$ imes$ 支持

 优化方法,即用优化方法求解函数距离零点最近,具体方法为信赖域方法。

包含预处理共轭梯度(PCG)、狗腿(dogleg)算法和Levenberg-Marquardt算法等

fzero由初值得到一个实根由初值得到一个实根支持$ imes$ $ imes$ 支持

一维解非线性方程方法,二分法、二次反插和割线法的混合运用

具体原理见数值求解非线性方程的二分法、不动点迭代和牛顿法插值方法 

roots支持,得到全部数值解$ imes$支持$ imes$$ imes$ 支持

特征值方法,即将多项式转化友矩阵(companion matrix)

然后使用求矩阵特征值的算法求得所有解(那是另外一个问题了)

   也就是说,之前写了几篇关于非线性求解的,如二分法、牛顿法(参见二分法、不动点迭代、牛顿法)、二次反插法(参见插值法),只有一个库函数用了类似的方法。

各函数用法详解

1. (符号/数值)解方程(组)函数 solve

  官方参考页:Equations and systems solver - MATLAB solve

  solve 是基本的用于符号解方程的内置函数,返回类型为符号变量矩阵($m imes 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

  非线性方程(组):MATLAB内置函数 solve, vpasolve, fsolve, fzero, roots [MATLAB]第1张

  输出内容中,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);

  输出图片:

非线性方程(组):MATLAB内置函数 solve, vpasolve, fsolve, fzero, roots [MATLAB]第2张

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

  

免责声明:文章转载自《非线性方程(组):MATLAB内置函数 solve, vpasolve, fsolve, fzero, roots [MATLAB]》仅用于学习参考。如对内容有疑问,请及时联系本站处理。

上篇Linux cat查看文件,查找关键字(grep),统计(wc -l)SpringBoot源码深度解析下篇

宿迁高防,2C2G15M,22元/月;香港BGP,2C5G5M,25元/月 雨云优惠码:MjYwNzM=

相关文章

BUAA_2019_MATLAB基础与应用_期末复习纲要

Matlab复习提纲 一、概述 1. Matlab(Matrix Laboratory)概述 1980年,由美国的 Clever Moler 博士开发; 是一款 科学与工程计算软件; 第四代智能计算机语言。 2. 功能与特点 开放性强、可扩展性强,兼容性强,直观灵活; MATLAB提供了丰富的矩阵运算处理功能,是基于矩阵运算的处理工具; 矩阵运...

《量化投资:以MATLAB为工具》连载(1)基础篇-N分钟学会MATLAB(上)

http://blog.sina.com.cn/s/blog_4cf8aad30102uylf.html 《量化投资:以MATLAB为工具》连载(1)基础篇-N分钟学会MATLAB(上) 《量化投资:以MATLAB为工具》简介 《量化投资:以MATLAB为工具》是由电子工业出版社(PHEI)下属旗舰级子公司——北京博文视点资讯有限公司出版的《量化投资与对冲...

深度自编码器(Deep Autoencoder)MATLAB解读

深度自编码器(Deep Autoencoder)MATLAB解读作者:凯鲁嘎吉 - 博客园http://www.cnblogs.com/kailugaji/ 这篇文章主要讲解Hinton在2006年Science上提出的一篇文章“Reducing the dimensionality of data with neural networks”的主要思想与M...

matlab初探寻

1 matlab 2 <iframe src="http://player.bilibili.com/player.html?aid=74994893&cid=128293306&page=1" scrolling="no" border="0" frameborder="no" framespacing="0" allowfu...

matlab调教日记 --- debug篇

节运行 在进行matlab测试的时候,要尽量采用节运行的功能,避免其他代码对调试部分代码进行的干扰。 节运行方式如下图,在注释中加入%%既可以进行分节 点击matlab 编辑器工具栏中的运行节 即可仅运行该节。 若不喜欢该方法,也可以使用将其他代码全部注释,或者使用matlab新推出的交互脚本方式进行调试,但交互脚本方式似乎运行速度较为缓慢,适合小型...

Matlab/Simulink仿真中如何将Scope转化为Figure?

1.只需要在运行仿真后,在命令窗口内输入: set(0,'ShowHiddenHandle','on'); set(gcf,'menubar','figure'); scope最上方会出现一个菜单栏,选择Tools->Edit Plot,即可修改图像所有属性。 2.双击Scope->Parameters->Data History...