【原】14.1 解方程 (数值方法)

数值方法的介绍

计算机在科学研究中的一个重要应用是找到那些没有解析解（即不能用多项式和标准数学函数写出来的解）的数学问题的数值解。在本章中，我们简要介绍了一些应用了高度发展的数值方法的领域，例如解非线性方程、求解积分和解微分方程。

14.1 解方程

在本节中，我们将讨论如何数值地解一元方程。通常，可以通过表达式 来描述这个问题，也就是说，我们希望找到方程的根（或者说解）。这个过程也可以描述为寻找 的零点。对于任意的 ，没有一般的方法可以解析地找到它的根。

• 14.1.1牛顿法

• 二分法

• 复根

• 14.1.2 二分法

• 14.1.3 fzero

• 14.1.4 roots

• solve

• 补充

14.1.1牛顿法

牛顿法（Newton's method）可能是最容易实现的数值方法之一，用于求解方程，在前几章节中已经简要介绍过。它是一种迭代方法，意味着它会反复尝试改进对根的估计。如果 是对根的近似值，我们可以使用图中的直角三角形来将其与下一个近似值 联系起来。

其中， 。解得

实施牛顿法的结构计划如下：

1. 输入起始值 和所需的相对误差 。
2. 当相对误差 时，重复以下步骤直到达到 （可以自定义迭代次数）：
• 计算。
• 打印 和 。
3. 停止。

在第二步中限制迭代次数是必要的，因为该过程有可能不收敛。

尝试使用牛顿法找到方程 的非零根。这个例子可能会有些困难。如果你遇到困难，那么你就会发现牛顿法的一个严重问题：它只有在起始猜测值足够接近根时才会收敛。由于"足够接近"取决于 的性质和根的性质，所以显然可能会遇到困难。唯一的解决办法是对初始猜测值进行一些有智慧的试验和误差调整工作，最好的方法是通过绘制 的草图或使用MATLAB绘图，从而更容易做出准确的初始猜测。

二分法

如果牛顿法无法找到根，可以使用下面讨论的二分法（Bisection method）。详细步骤如下：

1. 选择一个区间，使得 和 的符号不同。这个区间应该包含待求解的根。

2. 将区间等分为两个子区间：和，其中 为子区间的中点。

3. 计算 。

4. 检查 是否等于零。如果是，那么 就是根，算法结束。

5. 如果 和 具有相同的符号，将子区间替换为新的区间，并返回第2步。

6. 如果 和 具有相同的符号，将子区间替换为新的区间，并返回第2步。

7. 重复步骤2到步骤6，直到找到满足精度要求的近似根。

二分法的优点是可以保证收敛到根，并且在每次迭代中都减少根的搜索区间。然而，相对于牛顿法，它的收敛速度较慢。

---

复根

牛顿法也可以找到复数根，但前提是起始猜测值必须是复数。

对于方程，假设以复数 作为初始值 ，使用以下脚本可以找到一个复数根。将disp( [x f(x)] )替换为disp( x )以输出根的值：

steps = 0; % 迭代计数器
x = input('初始猜测值：'); % 根的估计值
re = 1e-8; % 所需相对误差
myrel = 1;
while myrel > re & (steps < 20)
    xold = x;
    x = x - f(x)/df(x);
    steps = steps + 1;
    disp(x)
    myrel = abs((x-xold)/x);
end
if myrel <= re
    disp('找到零点')
    disp(x)
else
    disp('未找到零点')
end

运行该脚本，输出的结果应为：

0.0769 + 0.6154i
-0.5156 + 0.6320i
-0.4932 + 0.9090i
-0.4997 + 0.8670i
-0.5000 + 0.8660i
-0.5000 + 0.8660i
找到零点
-0.5000 + 0.8660i

由于复数根以复共轭对的形式出现，另一个根是。

14.1.2 二分法

再次考虑解方程 的问题，其中 。

在双分法中，我们尝试通过观察或试错找出两个值 和 ，使得 和 具有不同的符号，即。如果我们可以找到这样的两个值，根必定位于它们之间的区间内，因为 在这个区间内改变了符号。在这个例子中，和是可以的，因为 ，并且 。在双分法中，我们通过 来估计根，其中 是区间的中点，即：

然后，如果与具有相同的符号（如图中所示），根明显位于 和   之间。因此，我们需要重新定义区间的左端点，并将其值设为 ，即新的 值为 。否则，如果 和 有不同的符号，我们将新的   值设为 ，因为根必定位于 和 之间。重新定义 或   后，根据公式（14.1）再次对新的区间进行双分，并重复这个过程，直到 和   之间的距离足够小为止。

二分法的一个很好的特点是，在开始之前，我们可以计算出需要多少次二分法才能达到一定的精度，给定 和 的初始值。假设我们以 和 开始。经过第一次二分， 的最大可能误差 为 ，因为我们估计根位于区间的中点。最坏的情况是根实际上位于 或 ，此时误差是 。继续下去，经过 次二分，最大可能误差 由 给出。如果我们希望确保这个误差小于某个指定的误差E，我们必须确保 满足不等式 ，即，

为了解出 ，我们可以对不等式进行求解，得到

根据不等式（14.2）， 是二分法进行的双分次数，必须是一个整数。最小的整数 ，使其超过不等式（14.2）的右边，可以作为保证给定精度 所需的最大二分次数。

以下方案可用于编写二分法的程序。它适用于任何在两个值 和 之间改变符号（向任一方向）的函数 ，这两个值必须事先由用户确定。

1. 输入 和 。
2. 初始化 和 。
3. 根据不等式（14.2）计算最大二分次数
4. 重复 次：
• 根据公式（14.1）计算 。
• 如果 ，则令 ，否则令 。
5. 显示根 。
6. 停止。

我们假设这个过程不会精确地找到根;对于实变量，这种情况发生的几率是无穷小的。

二分法的主要优点是，如果可以找到两个起始值 ，使得函数在它们之间改变符号，那么它可以保证找到一个根。你还可以预先计算出实现给定精度所需的二分次数。与牛顿法相比，它效率较低。连续的二分不一定会更接近根，而这通常在牛顿法中发生。实际上，有趣的是，在同一个函数上比较两种方法，看看二分法比牛顿法需要多少步骤。例如，对于解方程，二分法需要21步才能达到牛顿法5步达到的相同精度。

14.1.3 fzero

fzero是MATLAB中的一个函数，用于数值求解非线性方程的根。它使用的是迭代方法，类似于牛顿法。

使用fzero函数，你需要提供一个函数句柄来表示你要求解的方程，以及一个初始猜测值。fzero将通过迭代方法，尝试找到方程的根。

以下是使用fzero函数的基本语法：

x = fzero(fun, x0)

• fun是指向你要求解的方程的函数句柄。这个函数应该接受一个变量（通常是x）作为输入，并返回方程的值。
• x0是方程的初始猜测值。

fzero函数将尝试寻找方程在初始猜测值周围的根。如果成功找到根，它将返回一个近似的根值。如果无法找到根，或者找到的根不稳定，则可能会返回一个警告或错误。

示例一：

fun = @(x) x^2 + x + 1;
x0 = 1; % 初始猜测值
x = fzero(fun, x0);
disp(x); % 显示根的值

以上示例中，我们定义了一个匿名函数fun，表示方程。我们将初始猜测值设置为1，并使用fzero函数求解方程的根。最后，我们将根的值显示出来。

请注意，fzero函数只能找到一个根。如果方程有多个根，你需要提供不同的初始猜测值，或者使用其他的数值求解方法。此外，fzero函数只能用于实数根的求解，不能直接用于复数根的求解。如果要求解复数根，你需要使用其他的方法或工具。

---

示例二 要找到函数 的零点，你可以使用MATLAB中的fzero函数。然而，fzero仅设计用于寻找实根，不能直接找到复数根。

fun = @(x) x^3 + x - 3;
a = 1; % 初始猜测值
x = fzero(fun, a);
disp(x); % 显示找到的零点

如果你想要找到函数的复数根，你需要使用其他方法或工具。MATLAB提供了多种寻找复数根的数值方法，例如roots函数或使用符号数学中的solve函数。这些方法可以处理复数根和实根。

14.1.4 roots

MATLAB中的函数文件roots(c)用于找到系数为向量c的多项式的所有根（零点）。

要使用roots函数，你需要提供一个长度为n+1的向量c，其中n是多项式的次数。向量c的元素表示多项式的系数，按降序排列。

以下是使用roots函数的基本语法：

r = roots(c)

• c是多项式的系数向量。
• r是包含多项式的所有根的向量。

roots函数将计算多项式的所有根，并返回一个向量，其中包含这些根。返回的根可以是实数或复数，取决于多项式的特性。

下面是使用roots函数找到多项式根的示例：

c = [1, 0, 1, -3]; % 多项式系数向量为 [1, 0, 1, -3]
r = roots(c);
disp(r); % 显示多项式的根

  -0.6067 + 1.4506i
  -0.6067 - 1.4506i
   1.2134 + 0.0000i

在这个示例中，我们定义了一个多项式的系数向量 c，其中多项式为 。然后，我们使用roots函数找到这个多项式的所有根，并将根的向量显示出来。

需要注意的是，根的顺序在返回的向量中是任意的，并且可能包含重复的根。如果多项式有重复的根，它们可能以不同的多项式顺序出现多次。

另外，对于高次数的多项式，求解所有的根可能会很耗时，甚至可能导致数值不稳定。在某些情况下，使用符号计算方法（如solve函数）可能更适合找到多项式的根。

---

solve

solve是MATLAB中用于符号计算的函数。它用于求解包含符号变量的方程或不定式。使用solve函数，你可以找到方程或不定式的解析解，并获得更详细的结果。

以下是使用solve函数的基本语法：

sol = solve(equations, variables)

• equations是一个方程或方程组，可以表示为符号表达式或字符向量的形式。每个方程都必须等于零（或与零相等）。
• variables是待解的符号变量，可以是一个符号向量或以逗号分隔的符号变量列表。

solve函数将尝试求解方程或方程组，并返回一个结构体数组sol，其中包含解析解。每个解都表示为一个结构体，包含各个变量的解析表达式。如果无法找到解，sol将为空。

以下是使用solve函数求解方程的示例：

syms x
eqn = x^3 + x - 3 == 0;
sol = solve(eqn, x);
disp(sol); % 显示方程的解析解

在这个示例中，我们首先定义了一个符号变量x。然后，我们定义了方程eqn，它表示为x^3 + x - 3 = 0。接下来，我们使用solve函数解决这个方程，并将解析解存储在sol中。最后，我们显示解析解。

运行以上代码，将得到方程的解析解。解析解是以符号形式表示的，可以在数学计算中使用或进一步处理。

solve函数还支持求解方程组、不定方程和其他符号计算问题。你可以参考MATLAB的文档以获取更多关于solve函数的详细信息和使用示例。

---

补充

以下是solve函数的常用的用法：

1. 求解单个方程

使用solve函数可以求解单个方程的解析解。以下是一个求解单个方程的示例：

syms x
eqn = x^2 - 4 == 0;
sol = solve(eqn, x);
disp(sol); % 显示方程的解析解

在这个示例中，我们首先定义了一个符号变量x。然后，我们定义了方程eqn，它表示为x^2 - 4 = 0。接下来，我们使用solve函数来求解这个方程，并将解析解保存在sol中。最后，我们显示解析解。

2. 求解方程组

solve函数还可以用于求解方程组，即多个方程组成的系统。以下是一个求解方程组的示例：

syms x y
eqn1 = x + y == 5;
eqn2 = x - y == 1;
sol = solve([eqn1, eqn2], [x, y]);
disp(sol); % 显示方程组的解析解

x: 3
y: 2

在这个示例中，我们定义了两个未知数x和y的方程组。方程eqn1表示为x + y = 5，方程eqn2表示为x - y = 1。然后，我们使用solve函数来求解这个方程组，并将解析解保存在sol中。最后，我们显示解析解。

注意，在方程组求解时，我们需要将方程组表示为一个符号表达式的向量，并且将未知数也表示为一个符号变量的向量。

3. 求解不定方程

除了求解方程和方程组，solve函数还可以用于求解不定方程。不定方程是指在方程中存在未知数，但并没有具体的等号关系。以下是一个求解不定方程的示例：

syms x
eqn = x^2 + 2;
sol = solve(eqn, x);
disp(sol); % 显示方程的解析解

-2^(1/2)*1i
2^(1/2)*1i

在这个示例中，我们定义了一个未知数x的不定方程。方程eqn表示为x^2 + 2。然后，我们使用solve函数来求解这个不定方程，并将解析解保存在sol中。最后，我们显示解析解。

注意，在不定方程求解时，方程中的未知数不需要等于零，而只需要找到满足方程的值。