MATLAB 求解线性和非线性方程

方程求解是科学与工程计算的基石。MATLAB 提供了世界一流的工具来求解各种类型的方程，从简单的线性系统到复杂的非线性方程组。其核心方法分为两大类：针对线性方程的直接法和针对非线性方程的迭代法。

• 线性方程 (Linear Equations): 形如 \(Ax = b\) 的方程。其特点是变量的次数均为一。如果解存在，通常可以通过直接的矩阵运算求得唯一解或通解。
• 非线性方程 (Nonlinear Equations): 包含变量的非线性项（如 \(x^2\), \(\sin(x)\), \(e^x\) 等）的方程。这类方程通常没有通用的解析解，必须使用迭代算法从一个初始猜测值开始，逐步逼近真实解。

线性方程求解

\ (mldivide) - 矩阵左除求解Ax=b

• 功能 这是 MATLAB 中求解线性方程组 \(Ax = b\) 的首选方法。它封装了高效、稳定的数值算法。

• 数学原理: 矩阵分解 从理论上讲，如果 \(A\) 是一个可逆的方阵，解可以表示为 \(x = A^{-1}b\)。然而，在数值计算中，直接计算矩阵的逆 \(A^{-1}\) 是一种非常低效且数值不稳定的方法。 取而代之的是，\ 运算符使用了更稳定、更快速的矩阵分解方法。它会自动检测矩阵 A 的性质并选择最优算法：

  • 三角矩阵: 使用前向或后向替换法，速度极快。
  • 对称/厄米正定矩阵: 使用乔列斯基 (Cholesky) 分解 （\(A = R^T R \,\)）。
  • 通用方阵: 使用LU 分解 （\(A = LU \,\)）。
  • 非方阵 (超定或欠定系统): 使用QR 分解 （\(A = QR \,\)） 来求解最小二乘解。

• 语法

  x = A \ b
  

• 语法说明

  • A: \(m \times n\) 的系数矩阵。
  • b: \(m \times 1\) 的结果列向量。
  • x: \(n \times 1\) 的解向量。

• 示例

  % 求解以下线性方程组:
  % 2x + y - z = 8
  % -3x - y + 2z = -11
  % -2x + y + 2z = -3
  
  % 1. 定义系数矩阵 A 和结果向量 b
  A = [ 2,  1, -1;
       -3, -1,  2;
       -2,  1,  2];
  b = [ 8; -11; -3];
  
  % 2. 使用反斜杠运算符求解
  x = A \ b;
  
  disp('解向量 x 是:');
  disp(x); % 结果应为 [2; 3; -1]
  
  % 3. 验证结果
  residual = A*x - b;
  disp('残差 (A*x - b):');
  disp(residual); % 结果应接近于零向量
  

linsolve - 通用线性系统求解器

• 功能 一个功能更丰富的线性方程求解器，可以利用矩阵的特殊结构（如上三角、对称正定）来进一步提高求解效率。

• 语法

  x = linsolve(A, b)
  x = linsolve(A, b, opts)
  

• 语法说明

  • linsolve(A, b): 与 A \ b 的基本功能相同。
  • opts: 一个结构体，用于告知求解器矩阵 A 的特殊属性，例如 opts.UT = true 表示 A 是上三角矩阵，linsolve 此时会直接使用后向替换法，避免了矩阵分解的开销。

• 示例

  A = [ 2,  1, -1; 0, -1,  2; 0,  0,  2]; % 这是一个上三角矩阵
  b = [ 8; -11; -3];
  
  % 告知 linsolve A 是上三角矩阵
  opts.UT = true;
  x = linsolve(A, b, opts);
  
  disp('使用 linsolve 和上三角选项的解是:');
  disp(x);
  

非线性方程求解

fzero - 单变量非线性函数求根

• 功能 寻找单个变量的非线性函数 \(f(x)\) 的根，即求解 \(f(x) = 0\)。

• 数学原理 fzero 结合了二分法 (Bisection Method)、割线法 (Secant Method) 和逆二次插值法 (Inverse Quadratic Interpolation)。

  • 如果提供一个包含根的区间 [x1, x2] (即 f(x1) 和 f(x2) 异号)，算法会保证找到根。它首先使用较慢但可靠的二分法来缩小范围，然后在根附近切换到收敛更快的割线法或逆二次插值法。
  • 如果只提供一个初始点 x0，算法会尝试从此点出发，向外搜索一个包含根的区间，然后再应用上述混合算法。

• 语法

  x = fzero(fun, x0)
  x = fzero(fun, [x1, x2])
  

• 语法说明

  • fun: 一个函数句柄，形如 @(x) ...。
  • x0: 一个初始猜测点。
  • [x1, x2]: 一个包含根的区间，要求 fun(x1) 和 fun(x2) 的符号必须相反。

• 示例

  % 求解方程 cos(x) = x
  % 首先，将其重写为 f(x) = cos(x) - x = 0
  
  % 1. 定义函数句柄
  my_fun = @(x) cos(x) - x;
  
  % 2. 提供一个初始猜测值
  x0 = 0.5;
  
  % 3. 使用 fzero 求解
  root = fzero(my_fun, x0);
  
  disp(['方程 cos(x) = x 的解是: x = ', num2str(root)]); % 结果约为 0.7391
  

fsolve - 多变量非线性方程组求解

• 功能 求解多变量的非线性方程组 \(F(x) = 0\) (需要 Optimization Toolbox)。

• 数学原理 fsolve 主要使用信赖域 (Trust-Region) 算法，特别是信赖域狗腿算法 (Trust-Region Dogleg)。其核心思想是：

  1. 从一个初始猜测值 \(\bm{x}_0\) 开始。
  2. 在每一步迭代中，在当前点 \(\bm{x}_k\) 的一个“信赖域”（通常是一个球形区域）内，构建一个简化的模型（通常是线性或二次模型）来近似真实的非线性函数 \(F(\bm{x})\)。
  3. 在信赖域内求解这个简化模型的根，得到一个“试探步长” \(\bm{s}_k\)。
  4. 评估这个试探步长的效果。如果它使得实际的残差 \(||F(\bm{x}_k + \bm{s}_k)||^2\) 显著下降，则接受这一步，并可能扩大信赖域；否则，缩小信赖域，重新求解。
  5. 重复此过程，直到解收敛。

• 语法

  x = fsolve(fun, x0)
  x = fsolve(fun, x0, options)
  

• 语法说明

  • fun: 一个函数句柄，它接受一个向量 x 并返回一个向量 F(x)，F 是方程组的残差。
  • x0: 初始猜测值的向量。
  • options: 由 optimoptions 创建的选项结构体，用于控制算法细节（如算法选择、显示迭代过程等）。

• 示例

  % 求解方程组:
  % 2x - y = exp(-x)
  % -x + 2y = exp(-y)
  
  % 1. 定义一个函数，输入为 x=[x(1); x(2)]，输出为 F=[F(1); F(2)]
  fun = @(x) [2*x(1) - x(2) - exp(-x(1));
               -x(1) + 2*x(2) - exp(-x(2))];
  
  % 2. 提供初始猜测值向量
  x0 = [-5; -5];
  
  % 3. 调用 fsolve
  options = optimoptions('fsolve', 'Display', 'iter'); % 显示迭代过程
  solution = fsolve(fun, x0, options);
  
  disp('方程组的解 [x, y] 是:');
  disp(solution'); % 结果约为 [0.5671, 0.5671]
  

符号方程求解 (Symbolic Math Toolbox)

solve - 通用符号方程求解器

• 功能 尝试寻找代数方程或方程组的精确解析解。

• 数学原理 solve 利用一个庞大的代数规则和算法库来求解方程。

  • 多项式方程: 使用求根公式，如二次方程的

    \[x = \frac{-b \pm \sqrt{b^2-4ac}}{2a}\]

    或更高级的代数方法。

  • 线性方程组: 使用高斯消元法等。

  • 其他方程: 尝试使用各种代数代换、函数反演等技巧来分离变量。 如果找不到解析解，它会返回一个空的符号对象，或者在某些情况下会给出警告并尝试调用数值求解器 vpasolve。

• 语法

  sol = solve(eqn, var)
  sol = solve(eqns, vars)
  S = solve(eqns, vars, 'ReturnConditions', true)
  

• 语法说明

  • eqn, eqns: 一个或多个符号方程。使用 == 来定义方程。
  • var, vars: 一个或多个要求解的符号变量。
  • 'ReturnConditions', true: 返回解成立的条件以及参数的限制。

• 示例 (单个方程)

  syms x a b c
  % 求解二次方程 ax^2 + bx + c = 0
  eqn = a*x^2 + b*x + c == 0;
  sol_x = solve(eqn, x);
  
  disp('二次方程的通解是:');
  pretty(sol_x); % 以漂亮的格式显示公式
  

• 示例 (方程组)

  syms u v
  % 求解一个非线性方程组
  eqns = [u^2 - v^2 == 1, u + v == 2];
  vars = [u, v];
  
  sol = solve(eqns, vars);
  
  disp('方程组的解是:');
  disp(['u = ', char(sol.u)]); % u = 5/4
  disp(['v = ', char(sol.v)]); % v = 3/4