数值计算笔记之非线性方程的求解（二）迭代法

0、基本原理：逐次逼近

$f(x) = 0 \Leftrightarrow x = g(x)$

给定某初始解 $x_{0}$ ，按 $x_{k+1} = g(x_{k})$ $k = 0,1,2,\cdots$

若{ ${x_{k}}$ }收敛于 $x^{*}$ ，且 $g(x)$ 连续，则 $x^{*}$ 即为 $f(x) = 0$ 的解（ $\because \lim_{k\rightarrow \infty }x_{k+1} = g(\lim_{k\rightarrow \infty }x_{k})$ $\therefore x^{*} = g(x^{*})$ ）

迭代法收敛的充分条件

定理：若迭代函数 $g(x)$ 满足

在区间 $[a,b]$ 上 $g(x)^{'}$ 存在，且存在 $0<L<1$ ，使 $|g(x)^{'}|\leq L$
$\forall x\in [a,b],g(x)\in [a,b]$

则（1）、任取初值 $x_{0}\in [a,b]$ ，迭代法都收敛于方程的唯一实根。

（2）、 $|x^{*}-x_{k}\leq \frac{1}{1-L}|x_{k+1}-x_{k}|$ 用于误差估计

（3）、 $|x^{*}-x_{k}\leq \frac{L^{k}}{1-L}|x_{1}-x_{0}|$ 用于误差估计

例：用迭代法求方程 $2^{x}-4x = 0$ ，再 $[0,0.5]$ 内的根，要求 $|x_{k+1}-x_{k}| \leq \frac{1}{2} \times 10^{4}$ ，保留5位有效数字。

解：把原方程化为： $x = \frac{1}{4}\cdot 2^{x} = g(x)$

$g(x)^{'} = \frac{1}{4}\cdot ln2\cdot 2^{x}$ 存在，且 $|g(x)^{'}\leq \frac{1}{4}\cdot ln2\cdot 2^{0.5} < 1$
$g(0) \leq g(x)\leq g(0.5) \Rightarrow 0.25\leq g(x)\leq 0.353 \therefore g(x)\in [0,0.5]$

故 $g(x)$ 收敛。

n	Xn
0	0
1	0.25000
2	0.29730
3	0.30721
4	0.30933
5	0.30937

$\therefore x_{5} = 0.30937$ 为近似解。

一、牛顿迭代法解方程

原理：求解 $f(x)=0$ ，在函数f（x）图上任何一点画一条切线，然后确定切线与x轴的交点。这个方法需要一个初始值xo，此后的迭代公式为 $x_{n+1} = x_{n} - \frac{f(x_{n})}{f^{'}(x_{n})}$ 。具体理解：

如图，设 $x_{0}$ 是 $f(x)=0$ 的一个近似根，则由泰勒展开有： $f(x) \approx f(x_{0}) + f^{'}(x_{0})(x-x_{0}) = 0$ ，若 $f^{'}(x_{0}) = 0$ ，则可得 $x_{1}=x_{0}-\frac{f(x_{0})}{f^{'}(x_{0})}$ ，重复该过程，有 $x_{n+1} = x_{n} - \frac{f(x_{n})}{f^{'}(x_{n})}$ 。

matlab代码：

clear;
clc;
% fun = 2* x^3 - 4* x^2 + 3*x - 6;
x0 = 0;
error = 1;
f = subs(myfunction,x0);
df = subs(diff(myfunction),x0);
while error > 0.0001f = double(subs(myfunction,x0));          %加上double是为了显示结果为小数df = double(subs(diff(myfunction),x0));
%     f = sin(x0)+ x0 + 1;
%     df = cos(x0) + 1;root = x0 - f/df ;         %牛顿迭代公式error = abs(x0 - root);    %更新误差x0 = root;                 %更新初始值
end
format long
root

调用的 $myfunction$ 函数为：

function y = myfunction
syms x;
y = sin(x) + x + 1;
end

syms是定义符号变量的函数，求导需要用到，没有申明的话，无法求导。
matlab中subs()是符号计算函数，表示将符号表达式中的某些符号变量替换为指定的新的变量，常用调用方式为：

subs(S,OLD,NEW) 表示将符号表达式S中的符号变量OLD替换为新的值NEW。
$f = double(subs(myfunction,x0))$ ，加上double是为了显示结果为小数

C++代码：

/* 牛顿迭代法 求方程根 */
#include <iostream>
#include"function.h"
using namespace std;int main()
{double x0, err, e, root;cout << "请输入初始值x0, x0 = "  ;cin >> x0;cout << "请输入误差限e, e = ";cin >> e;err = 10;                    //定义一个较大的数，使其进入循环while (err > e) {int n = 1;                   //迭代次数double x = x0;/*   double f = 2.0 * pow(x, 3) - 4 * pow(x, 2) + 3 * x - 6;double df = 6 * pow(x, 2) - 8 * x + 3;*/double f = Original_equation(x);double df = Integral_equation(x);x = x0 - f / df;root = x;err = fabs(x0 - root);x0 = x;n++;if (n > 150) {cout << "迭代次数太多，请更换算法，已推出程序！" << endl;break;}}cout << "方程的根为：root = " << root << endl;//cout << "Hello World!" << endl;
}

$function.h$ 头文件为：

#pragma once
double Original_equation(double &x) {double y = 2.0 * pow(x, 3) - 4 * pow(x, 2) + 3 * x - 6;return y;
}double Integral_equation(double &x) {double dy = 6 * pow(x, 2) - 8 * x + 3;return dy;
}

结果为：

但是，作为计算一般性函数零解问题的算法，它有三个严重的缺陷。

要求函数 $f^{'}(x)$ 必须是光滑的；
可能遇到导数 $f^{'}(x)$ 不方便计算的困难；
初始解必须靠近准确解。

至于优点，就是对于计算平方根的问题，牛顿法特别简洁有效。

二、弦截法（弦割法）

有的时候，求导是件很麻烦的事，所以弦截法就出来了。

首先了解一下差商：称 $\frac{f(x_{k+1}) - f(x_{k})}{x_{k+1} - x_{k}}$ 为 $f(x)$ 在 $x_{k}$ 、 $x_{k+1}$ 两点的差商。其实在 $x_{k}\rightarrow x_{k+1}$ 时，就是导数的定义。

原理：弦截法用最近两次选代解构造出的有限差分近似，替代牛顿法中的求导数计算，不同于牛顿法在f（x）曲线上某一点画切线的做法，它通过两个点画一条割线，下一个迭代解就是割线与x轴的交点。
割线法的迭代需要两个初始值，xo和x1，后续的迭代解按下面的公式计算：

1、定端点弦截法

$s_{k} = \frac{f(x_{k}) - f(x_{0})}{x_{k} - x_{0}}$ 、 $x_{n+1} = x_{n} - \frac{f(x_{k})}{s_{k}}$ (牛顿法中的 $f^{'}(x_{k})$ 被割线的斜率 $s_{k}$ 所替代)。

2、变端点弦截法

$s_{k} = \frac{f(x_{k}) - f(x_{k-1})}{x_{k} - x_{k-1}}$ 、 $x_{k+1} = x_{k} - \frac{f(x_{k})}{s_{k}} = x_{k} - \frac{f(x_{k})}{f(x_{k}) - f(x_{k-1})}\cdot (x_{k}-x_{k-1})$ (牛顿法中的 $f^{'}(x_{k})$ 被割线的斜率 $s_{k}$ 所替代)。

matlab代码：

clear;
clc;
f = @(x) x^3 - x -1;
x0 = 1;
x1 = 1.5;
k = 0;
while abs(x1 - x0) > 0.0001s = (f(x1) - f(x0)) / (x1 - x0);root = x1 - f(x1) / s;x0 = x1;x1 = root;k = k + 1;if k > 150fprintf("循环次数太多，推出程序"); disp('循环次数太多，推出程序');root = NULL;break;end
end
% format long
root

结果为：

C++代码：

#include <iostream>
#include"function.h"
using namespace std;int main()
{/*弦截法*/double x0, x1, error;double x ;double root;int n = 0;cout << "请输入初始值x0, x0 = ";cin >> x0;cout << "请输入初始值x1, x1 = ";cin >> x1;cout << "请输入误差限error, error = ";cin >> error;while (fabs(x1 - x0) > error) {double s = equation2(x1, x0);root = x1 - equation1(x1) / s;x0 = x1;x1 = root;n++;if (n > 150) {cout << "迭代次数太多，请更换算法，已推出程序！" << endl;break;}}cout << "方程的根为：root = " << root << endl;//cout << "Hello World!" << endl;
}

头文件：

#pragma once
double equation1(double &x) {double y = pow(x, 3) - x - 1;return y;
}
double equation2(double &x1, double &x0) {double s = (equation1(x1) - equation1(x0)) / (x1 - x0);return s;
}

结果为：