MATLAB-求解方程

迷失缠绵情愫中

求解方程通常有两种方法，符号求解和数值求解。

1. solve(Symbolic Math Toolbox)
通常在不确定方程是否有符号解的时候，推荐先使用 solve 进行尝试，因为 solve 相比于数值求解来说，它不需要提供初值，并且一般情况下能够得到方程的所有解。对于一些简单的超越方程，solve 还能够自动调用数值计算系统给出一个数值解。

solve 的常见调用形式：

1. sol = solve(eq)
   
2. sol = solve(eq,var)
   
3. sol = solve(eq1,eq2,…,eqn)
   
4. sol = solve(eq1,eq2,…,eqn,var1,var2,…,varn)

复制代码

eq 为符号表达式，var 为指定的要求解的变量。如果不声明要求解的变量(第一和第三种形式)，则 MATLAB 自动按默认变量进行求解，默认变量可以由 symvar(eq) 确定。

1. %例：求解方程组:  x+y = 1,  x-11y = 5
   
2. 
   
3. syms x y      %声明符号变量
   
4. eq1 = x+y-1;
   
5. eq2 = x-11*y-5;
   
6. sol = solve(eq1,eq2,x,y);
   
7. x = sol.x
   
8. y = sol.y

复制代码

solve 求得的解通过结构体的形式赋值给 sol，然后再通过 x = sol.x 和 y = sol.y 分别赋值给 x 和 y 。 也可以直接使用：

1. [x,y] = solve(eq1,eq2,x,y)

复制代码

进行求解，但需要注意，等式左边接收参数时应当按字母表进行排序，否则 MATLAB 不会自动识别你的参数顺序，比如：

1. [x,y] = solve(eq1,eq2,x,y)
   
2. [y,x] = solve(eq1,eq2,x,y)

复制代码

solve 会把答案按字母表进行排序后进行赋值，x 解赋值给第一个参数，y 解赋值给第二个参数，对于第二种形式，实际上最终结果是变量 y 存储了 x 的解而变量 x 存储了 y 的解。 上述情况另一种解决方案是用下面的方法指定输出顺序：

1. syms x y
   
2. [y,x] = solve(x + y == 2,x - y == 1, [y, x])

复制代码

由于是符号求解，有时候得到的解是一大串式子（符号求解无精度损失，所以 MATLAB 不会自动将答案转化为浮点数），这时候可以用 vpa 或者 double 函数将结果转换为单一的数，但需要注意的是 double 的结果为浮点数，vpa 的结果仍然是符号类型（即 sym 类型）。

另外，很多人习惯对于 solve 的参数采用字符型输入，这种方式有几个弊端。
首先就是程序的调试，一旦式子输入有误(最常见的就是括号的匹配)，则调试起来会非常困难，例如：

1. solve('10^(-4.74)*0.965*y/60000x/(10^(-4.74)+x)+0.1/36500+10^(-14)/x-x=0','10^(-3.2)*x+0.333/3000+8*10^((-3.2)*0.1+0.1/333*y','x','y')

复制代码

这时要去寻找式子输入错误会是一件很麻烦的事，MATLAB 也不会报告具体出错的地方。如果采用符号变量输入：

1. syms x y
   
2. eq1=10^(-4.74)*0.965*y/60000x/(10^(-4.74)+x)+0.1/36500+10^(-14)/x-x;
   
3. eq2=10^(-3.2)*x+0.333/3000+8*10^((-3.2)*0.1+0.1/333*y;
   
4. sol=solve(eq1,eq2,x,y)

复制代码

会给程序的调试带来许多便利。对于某些错误，MATLAB 会给出错误代码颜色的高亮， 命令行还能返回具体的错误信息。并且采用字符型输入时，变量的赋值不能传入方程：

1. %例：x+y*sin(x) = 1
   
2. 
   
3. y = 1;
   
4. sol = solve('x+y*sin(x)=1','x')

复制代码

MATLAB 会返回一个空解，而 sym 型输入：

1. syms x
   
2. y = 1;
   
3. eq = x+y*sin(x)-1;
   
4. sol = solve(eq,x)

复制代码

能够得到 sol = 0.51097342938856910952001397114508，其中的区别就在于 char 型输入尽管在 solve 前对 y 有一个赋值，但 solve 求解时依然会将 y 当作一个未赋值的常数。
最后，在今后的高版本 solve 将不支持 char 型参数输入，因此应该尽量放弃使用这种方法。



2. fzero
很多情况下， solve 并不能求得方程的解析解，这时可以采用数值法求解。
数值求解法包括 fzero 和 fsolve，其区别在于 fzero 只适用于求解一元函数零点，而 fsolve 适用于求解多元函数零点(包括一元函数)。
当求解一元函数零点时，推荐优先使用 fzero，原因是 fzero 求解一元方程往往更容易，它不仅支持提供初值的搜索，还支持在一个区间上进行搜索。

fzero 的常用形式：

1. x = fzero(fun,x0)
   
2. [x,fval] = fzero(fun,x0)

复制代码

其中 fun 为函数句柄， x0 为搜索初值， fval 为求解误差。

1. %例：一元方程 sin(x)+cos(x)^2 = 0
   
2. 
   
3. y = @(x)sin(x)+cos(x).^2;   
   
4. [x,fval] = fzero(y,1)  %1为搜索初值

复制代码

若方程有多个零点，fzero 只能根据你提供的初值求得最靠近初值的一个零点，如果希望求得多个零点，只能够通过改变初值来得到不同的零点。
对于初值的选取，目前来说没有什么比较好的办法，只能够通过分析方程的性质，或者通过作图的方法去寻找一个比较靠近零点的初值。另外，fzero 能够提供区间搜索，注意区间两端的端点函数值符号需要反向：

1. y = @(x)sin(x)+cos(x).^2;  
   
2. [x,fval] = fzero(y,[-1 1])     %fzero在[-1,1]这个区间进行搜索

复制代码

建议尽量用区间搜索的方式来求解，因为这种方法比单纯的提供一个初始值的运算速度要快一些。而且新版本的 MATLAB 中关于此函数还有多个参数的形式。



3. fsolve(Optimization Toolbox)
fsolve 可以求解多元方程，用法和 fzero 类似。 fsolve 的常用形式：

1. x = fsolve(fun,x0)
   
2. [x,fval] = fsolve(fun,x0)

复制代码

其中 fun 为函数句柄， x0 为搜索初值， fval 为求解误差。

1. %例：求解方程组 x+y = 1, x-11y = 5
   
2. 
   
3. eq = @(x)[x(1)+x(2)-1;x(1)-11*x(2)-5];
   
4. [sol,fval] = fsolve(eq,[1,1])

复制代码

这里对于方程的的输入需要采用矩阵的形式，其中 x(1) 代表 x ， x(2) 代表 y 。有时候变量较多时可能会容易混淆，这里提供另一种方法，用符号变量表达方程再利用 matlabFunction 转化为函数句柄：

1. syms x y
   
2. eq1 = x+y-1;
   
3. eq2 = x-11*y-5;
   
4. eq1 = matlabFunction(eq1);   %将符号函数转化为函数句柄
   
5. eq2 = matlabFunction(eq2);
   
6. eq = @(x)[eq1(x(1),x(2)); eq2(x(1),x(2))];
   
7. [sol,fval] = fsolve(eq,[1,1])

复制代码

结果与之前相同，但不容易出错。求得的解以矩阵形式返回给 sol ，即 sol 的第一个值求解的是 x(1) ，sol 的第二个值求解的是 x(2) 。
fsolve 要求求解的函数必须是连续的，而且成功求解时，fsolve 只能给出一组根。缺省情况下，trust-region dogleg 算法只能用于系统方程是方阵的情况，而 Levenberg-Marquardt 算法没有此限制。



4. vpasolve(Symbolic Math Toolbox)
最后再补充一个数值解法 vpasolve，vpasolve 是 R2012b 引进的函数，可以求解一元或多元函数零点。相比于 fzero 和 fsolve 来说，vpasolve 最大的一个优点就是不需要提供初值，并且能够自动搜索指定范围内的多个解。

vpasolve 调用形式：

1. S = vpasolve(eqn)
   
2. S = vpasolve(eqn,var)
   
3. S = vpasolve(eqn,var,init_guess)
   
4. ___ = vpasolve(___,Name,Value)

复制代码

其中 eqn 是符号方程，var 为待求解变量，也可以不提供(第一种调用形式，默认求解变量由 symvar(eqn) 求得)， init_guess 为搜索初值，Name，Value 为选项控制。

1. %例：对于多项式方程，vpasolve 能够给出所有解
   
2. 
   
3. syms x
   
4. vpasolve(4*x^4 + 3*x^3 + 2*x^2 + x + 5 == 0, x)
   
5. 
   
6. ans =
   
7. 
   
8. - 0.88011 - 0.76332i
   
9.    0.50511 + 0.81599i
   
10.    0.50511 - 0.81599i
    
11. - 0.88011 + 0.76332i

复制代码

对于非多项式方程，vpasolve 只能给出一个解：

1. syms x
   
2. vpasolve(sin(x^2) == 1/2, x)
   
3. 
   
4. ans =
   
5. 
   
6. -226.94

复制代码

这时可以提供搜索初值，来搜寻其它零点：

1. syms x
   
2. vpasolve(sin(x^2) == 1/2, x,100)
   
3. 
   
4. ans =
   
5. 
   
6. 99.996

复制代码

可以指定搜索范围，但不同于 solve，solve 指定求解范围是用 assume 函数，vpasolve 则是直接在输入参数中指定：

1. syms x
   
2. vpasolve(x^8 - x^2 == 3, x, [-Inf Inf]) %实数范围内求解

复制代码

最后，vpasolve 一个很强大的用法，将 ‘random’ 选项设置为 true 可以直接搜索指定范围内不同解：

1. syms x
   
2. f = x-tan(x);
   
3. for n = 1:3
   
4.   vpasolve(f,x,'random',true)
   
5. end

复制代码