方程数值解法--固定点迭代法与牛顿法原理与c++实现

发布时间:2024年01月04日

方程数值求根的常用方法有二分法、固定点迭代法、牛顿法等。这里介绍后面两种方法。

首先定义问题:对于函数 y = f ( x ) y=f(x) y=f(x),在求解零点(根)时,有:
f ( x ) = 0 (1) f(x)=0\tag{1} f(x)=0(1)

不动点迭代法:

我们可以将上式转化为:
x = g ( x ) (2) x=g(x)\tag{2} x=g(x)(2)
的形式。上式即为不动点迭代的方程。不动点迭代的具体做法如下:

(1)给定初值 x k = x 0 x_k=x_0 xk?=x0?

(2)令 x k + 1 = g ( x k ) x_{k+1}=g(x_k) xk+1?=g(xk?)

(3)如果不满足迭代条件继续执行(2),否则终止迭代,则最终方程的根 x ? = x k + 1 x*=x_{k+1} x?=xk+1?.

收敛性条件:

要求 g ( x ) g(x) g(x)的导数在(0,1)之间。证明过程可以看下面文章:

[1](数值分析学习笔记(二) - 知乎 (zhihu.com))

2

[3](非线性方程组求解方法合集——不动点迭代法(包含算法流程及代码) - 知乎 (zhihu.com))

例子: 2 x 3 ? x ? 1 = 0 2x^3-x-1=0 2x3?x?1=0

g ( x ) g(x) g(x)的形式是多样的,例如 x = 2 x 3 ? 1 x=2x^3-1 x=2x3?1,或者 x = x + 1 2 3 x=\sqrt[3]{\frac{x+1}{2}} x=32x+1? ?.对于第一种形式,显然不满足收敛性,所以选择第二种形式:

void fixedPointIter(const double init_value,const int max_iter, double& final_value)
{
	double x = init_value;
	int iter_count = 0;
	double x_1 = x;
	while (true)
	{
		
		//x_1 = 2.0 * std::pow(x,3)  - 1.0;
		x_1 = std::pow((x + 1) * 0.5, 1.0 / 3.0);
		
		std::cout << x_1 << "\n";

		if (iter_count>max_iter||std::abs(x_1-x)<1e-3)
		{
			break;
		}
		else
		{
			x = x_1;
			++iter_count;
		}
	}
	final_value = x_1;
}
int main()
{
	double x = 0.0;
	fixedPointIter(0.0, 50, x);
	std::cout<<x << "\n";
}

牛顿迭代法

对于式(1),理所当然可以写成: x = x + f ( x ) x=x+f(x) x=x+f(x),但是不满足导数小于1的收敛条件,因此,我们将其写为更好地形式,即令 g ( x ) = x ? f ( x ) f ′ ( x ) g(x)=x-\frac{f(x)}{f'(x)} g(x)=x?f(x)f(x)?,所以有牛顿法:
x k + 1 = x ? f ( x k ) f ′ ( x k ) (3) x_{k+1}=x-\frac{f(x_k)}{f'(x_k)}\tag{3} xk+1?=x?f(xk?)f(xk?)?(3)
更多理论说明[0](02非线性方程求根:牛顿法 - 知乎 (zhihu.com)),当然也可以从几何角度解释[1]((99+ 封私信 / 80 条消息) 如何通俗易懂地讲解牛顿迭代法求开方(数值分析)? - 知乎 (zhihu.com))

例子: 2 x 3 ? x ? 1 = 0 2x^3-x-1=0 2x3?x?1=0,code:

void NewtonIter(const double init_value, const int max_iter, double& final_value)
{
	double x = init_value;
	int iter_count = 0;
	double x_1 = x;
	
	while (true)
	{
		double f_x = 2 * x * x * x - x - 1;
		double f_x_derivative = 6 * x * x - 1;
		x_1 = x - f_x / (f_x_derivative );
		std::cout << "newton x:" << x_1 << "\n";

		if (iter_count > max_iter || std::abs(x_1 - x) < 1e-3)
		{
			break;
		}
		else
		{
			x = x_1;
			++iter_count;
		}
	}
	final_value = x_1;
}
int main()
{
	double x = 0.0;
	NewtonIter(0.0, 50, x);
	std::cout << "newton final x:" << x << "\n";
}

文章来源:https://blog.csdn.net/xinxiangwangzhi_/article/details/135392190
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。