数值分析——二分法和牛顿迭代（Bisection Method & Newton‘s Method）

本系列整理自博主21年秋季学期本科课程 数值分析I 的编程作业，内容相对基础，参考书: David Kincaid, Ward Cheney - Numerical Analysis Mathematics of Scientific Computing (2002, Americal Mathematical Society)

目录

适用对象

二分法

（1）算法描述

（2）代码

牛顿迭代法

（1）算法描述

 （2）代码

（3）定理

1. Theorem on Newton's Method

2. Theorem on Newton's Method for Convex Function

拓展

（1）弦方法

（2）函数迭代（Functional Iteration）

1. 压缩映射定理（Contractive Mapping Theorem）

2. 误差分析

总结

---

适用对象

求非线性方程的数值解，即寻找x使得目标函数f(x)=0

二分法

（1）算法描述

要求f(x)在区间[a,b]上连续，且f(a)f(b) < 0.

令c=a+(b-a)/2，检验f(a)f(c) < 0是否为真，若真，则令b=c；若否，且f(a)f(c)!=0，则令a=c.

重复上述过程，直到f(a)f(c)=0即f(c)=0时停止。数值上，一般要求|f(c)|的值小于某个数，或区间长度小于某个数时，停止程序。

该算法线性收敛于方程的根。

（2）代码

#include<iostream>
#include<math.h>
#include <iomanip>
using namespace std;
int M 50;//set the step M large enough
double f(double x){}//the target function

int sign(double x) {
	if (x > 0)
		return 1;
	else if (x == 0)
		return 0;
	else return -1;
}
int bisection(double a, double b) {
	//apply bisection on function f on interval [a,b]
	double left = a, right = b, length = mid = 0, valuel = f(left), valuer = f(right), valuem = 0;
	cout << "---------------------------------------------------------------------------" << endl;
	cout << setw(2) << " k" 
		 << setw(10) << "left"<< setw(13) << "f(left)" 
		 << setw(10) << "right" << setw(15) << "f(right)" 
		 << setw(10) << "mid" << setw(13) << "f(mid)" << endl;
	cout << "---------------------------------------------------------------------------" << endl;
	if (sign(valuel) == sign(valuer)) {
		length = right - left;
		mid = left + length / 2;
		valuem = f(mid); valuel = f(left); valuer = f(right);
		cout << setw(6) << 1 
			 << setw(10) << left<< setw(15) << valuel 
		     << setw(10) << right << setw(15) << valuer 
			 << setw(10) << mid << setw(15) << valuem<< endl;
		cout << "---------------------------------------------------------------------------" << endl;
		cout << "f(x) signs same at both ends of the interval ["< a < <", "< < b < <"], cannot use bisection algrithm to find the root." << endl;
		return -1;
	}
	for (int k = 1; k <= M; k++) {
		length = right - left;
		mid = left + length / 2;
		valuem = f(mid); valuel = f(left); valuer = f(right);
		cout << setw(6) << k
			 << setw(10) << left << setw(15) << valuel
			 << setw(10) << right << setw(15) << valuer
			 << setw(10) << mid << setw(15) << valuem << endl;
		if (abs(length) < 0.0001 || abs(valuem) < 0.0001)//stop when the value of function is smaller than 0.0001
			break;
		if (sign(valuem) != sign(valuel))
			right = mid;
		else
			left = mid;
	}
	cout << "---------------------------------------------------------------------------" << endl;
	if (abs(length) < 0.0001 && abs(valuem) > 100)//this determination condition should be modified according to different functions
		cout << "f(x) has no root on interval [" << a << ", " << b << "]，and it tend to infinity at the point x=" << mid  << endl;
	else
		cout << "the root of f(x) on interval [" << a << "," << b << "] is x = " << mid << endl;
	return 1;
}

牛顿迭代法

（1）算法描述

令r为f的零点，x趋向于r，由Taylor定理：0=f(r)=f(x+h)=f(x)+h*f'(x)+O(h^2)，其中O(h^2)可以忽略，推出h= - f(x)/f'(x)，从而零点 r=x+h=x - f(x)/f'(x)

牛顿迭代公式：

该算法二次收敛于方程的根。

 （2）代码

#include<iostream>
#include<math.h>
#include <iomanip>
using namespace std;
int M = 50;//the maximum of operation step
double f(double x){}//the target function
double df(double x){}//the derivative of f

double Newton_Method(double x) {//x is the initial value
	cout << " Newton's method" << endl;
	cout << "-----------------------------------------------" << endl;
	cout << setw(9) << " k" << setw(13) << "x_k" 
         << setw(15) << "f(x_k)" << setw(15) << "df(x_k)"
		 << resetiosflags(ios::left) << endl;
	cout << "-----------------------------------------------" << endl;
    double y = f(x);
	for (int k = 1; k <= M; k++) {
		cout << setw(7) << k << setw(13) << x 
             << setw(15) << y << setw(15) << df(x)<< endl;
		if (abs(y) == 0)
			break;
		x = x - y / df(x);
		y = f(x);
	}
	cout << "-----------------------------------------------" << endl;
	cout << "the root of the function f is x = " << x << endl;
	return 0;
}

（3）定理

1. Theorem on Newton's Method

对于牛顿迭代，若 f'' 连续，r 为 f 单根，则存在 r 的领域 B 和常数 c，当初值 x0∈B 时，xn 将稳定趋于 r，且二阶收敛

若 r 为 k 重根，将迭代公式修改为

且 xn 一阶收敛于 r

2. Theorem on Newton's Method for Convex Function

若 f 二次连续可导，单调，凸，且存在唯一零点，则任意初值，牛顿迭代都将收敛于该零点。

拓展

（1）弦方法

超线性收敛：

r 为零点，记 en=xn-r，由Taylor定理：f(xn)=f(r+en)=f(r)+en f'(r)+1/2 en^2 f''(r) +O(en^3)，推出 

f(xn)/en=f'(r)+1/2 en f"(r)+O(en^2)，同理f(xn-1)/en-1=f'(r)+1/2 en-1 f"(r)+O(en-1^2)

由(xn-xn-1)/(f(xn)-f(xn-1))≈1/f'(r)，

记C=1/2 f''(r)/f'(r)，则  

注：

弦方法收敛阶高于二分法，低于牛顿法。相比于牛顿迭代，弦方法不需要对函数求导，用“弦”来代替导数。

（2）函数迭代（Functional Iteration）

1. 压缩映射定理（Contractive Mapping Theorem）

设C为实数集R的闭子集，F：C->C压缩映射，则F有唯一不动点s，即任意初值x0∈C，xn收敛于s，且{xn}为Cauchy列。

2. 误差分析

记en=xn-s，若F'存在且连续，由中值定理，即

在不动点s处，设，且，则

对F(s+en)做Taylor展开，可得

即迭代法 q阶收敛于s

总结

算法不足：由于函数f使用全局定义，故程序一次只能操作一个函数，不能同时对多个函数求根。