非线性一元方程求解——弦截法、牛顿迭代法（C++）

一元方程

对于一元方程，如果要求f(x)=0的解，其过程大致包括如下三个问题：

根的存在性：是否有根，如果有，有几个？
根的分布：根分布区间；
求根的公式：如何从根的近似值逐步精确化？

求解方法

首先，对于某些特殊放入函数 $f (x)$ ，如果满足以下定理： $f (x)$ 在区间 $[a, b]$ 上连续， $f(a)\cdot f(b)<0$ ，则 $f (x) = 0$ 在区间 $[a, b]$ 上至少存在一个根。
然后就需要确定根的分布，此处有两种方法：

搜索法：将区间 $[a, b]$ 分割为n个子区间，子区间长度为 $\Delta x$ ，如果 $f(x_i)\cdot f(x_{i+1})<0(i=0,\cdots,n)$ ，则根在区间 $x_i,x_{i+1}]$ 内。
对分法：将区间 $[a, b]$ 从中点分开，分别计算 $f(a)\cdot f((a+b)/2)$ 和 $f(b)\cdot f((a+b)/2)$ ，判断哪个乘积小于0，继续对分对应区间，直至满足精度。

对于一般的情况，对分法的搜索速度较快，但是如果存在根与某个区间端点非常靠近的情况，搜索法反而速度更快。对分法还有一个问题，若区间内部存在多个根，可能会出现遗漏和算法报错等原因，因此需要根据实际需要选取对应的方法确定根的分布。
当初始区间较大时，对于上述两种方法其收敛速度都不太理想，因此引入下面的方法：

牛顿迭代法：将原方程进行泰勒展开，可以发现 $f(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{f''(\xi)}{2!}(x-x_0)^2$ ，因此有 $f(x)\approx f(x_0)+f'(x_0)(x-x_0)$ ，由此得到牛顿迭代式： $x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}$ 。可是牛顿迭代法的收敛性取决于 $x_0$ 的选取，因此使用牛顿迭代法之前还需要满足如下条件：
- $f(a)\cdot f(b)<0$ ;
- $f'(x)\neq 0$ ;
- $f^{''} (x)$ 存在且不变号;
- 取 $x_0\in [a,b]$ ，使得 $f''(x_0)\cdot f(x_0)> 0$ ;
弦截法：相比于牛顿迭代法，弦截法不需要使用导数，使用差商进行计算，如果 $f (x)$ 比较复杂，即使是近似计算导数也非常不方便。根据这个思路可以得到弦截法的迭代公式： $x_{k+1}=x_k-(x_k-x_{k-1})\frac{f(x_k)}{f(x_k)-f(x_{k-1})}$ 。弦截法的使用条件与牛顿迭代法相同，收敛速度也较牛顿法慢，并且需要两个初值。

代码实现

首先是牛顿法的实现：

#include<iostream>
#include<cmath>
#include <chrono>//计时模块
#include <iomanip>
#define DerivativeIntervalAccuracy 1e-6
//近似求导，传递函数进来
double f_derivative(double (*fx)(double),double x)
{return (fx(x + DerivativeIntervalAccuracy) - fx(x))/DerivativeIntervalAccuracy;
}//指定函数
double fx(double x)
{//计算一个复杂函数return pow(x,4)-11.7*log10(x+3)+exp(-x)-5;
}
double Newton_iterative(double (*fx)(double x),double x_0,double error)
{double x_1=0.0,x_2=x_0;//x_1表示x_k-1,x_2表示x_kdouble d_fx;int i = 0;do{x_1=x_2;d_fx=f_derivative(fx, x_1);x_2=x_1-fx(x_1)/d_fx;if(i++>10000){std::cout<<"迭代次数过多"<<std::endl;return 0.0;}//std::cout<<"x(k)="<<x_1<<" "<<"x(k+1)="<<x_2<<" "<<std::endl;}while(fabs(x_2-x_1)>error);//std::cout<<"迭代次数"<<i<<std::endl;return x_2;
}
int  main()
{double error=1e-6;double x;double x_0=100.0;//初始值double d_fx = 0.0;//导数auto start = std::chrono::steady_clock::now();x = Newton_iterative(fx, x_0, error);// 结束计时auto end = std::chrono::steady_clock::now();// 计算时间差auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start);// 输出时间差（以毫秒为单位）std::cout << "Elapsed time: " << duration.count() << " microseconds" << std::endl;std::cout<<"x="<<std::setprecision(10)<<x<<std::endl;return 0;
}

其次是弦截法的实现：

#include<iostream>
#include<cmath>
#include <iomanip>
#include <chrono>//计时模块
using namespace std;
//指定函数
double fx(double x)
{//计算一个复杂函数return pow(x,4)-11.7*log10(x+3)+exp(-x)-5;
}
double chord_truncation(double (*fx)(double x),double x_0,double x_1,double error)
{int i = 0;double x_2,x_3=x_0,x_4=x_1;do{x_2 = x_3;x_3 = x_4;x_4=x_3-(x_3-x_2)*(fx(x_3)/(fx(x_3)-fx(x_2)));if(i++>10000){std::cout<<"迭代次数过多"<<std::endl;return 0.0;}cout<<"x(k-1)="<<x_2<<" "<<"x(k)="<<x_3<<" "<<"x(k+1)="<<x_4<<endl;}while(fabs(x_3-x_2)>error);//std::cout<<"迭代次数"<<i<<std::endl;return x_4;
}
int  main()
{double error=1e-6;double x_1=10.0,x_2=0.0;//x_1表示x_k-1,x_2表示x_kauto start = std::chrono::steady_clock::now();double x=chord_truncation(fx,x_1,x_2,error);// 结束计时auto end = std::chrono::steady_clock::now();// 计算时间差auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start);// 输出时间差（以毫秒为单位）std::cout << "Elapsed time: " << duration.count() << " microseconds" << std::endl;cout<<"x="<<std::setprecision(15)<<x<<endl;return 0;
}