[数值分析]开方根算法 - 牛顿迭代法/二分法/割线法

// ConsoleApplication1.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include "math.h"
#include "stdlib.h"

#define err_limit	0.0001		//注意float的精度，不能要求太高，否则float无法区分
#define iter_limit	1000000

float NewtonSqrt(float c, int* iter_t);
float DivisionSqrt(float c, int* iter_t);
float SecantSqrt(float c, int* iter_t);

int _tmain(int argc, _TCHAR* argv[])
{
	float	num		= 123456;	//输入需要开方的数
	int		iter_t	= 0;		//迭代次数
	float	result	= 0;		//运算结果

	//牛顿迭代法求根
	result = NewtonSqrt(num, &iter_t);
	printf("牛顿迭代法求根\nNum:%f\nSqrtResult:%f\nTime:%d\n\n", num, result, iter_t);

	//二分法求根
	result = DivisionSqrt(num, &iter_t);
	printf("二分法求根\nNum:%f\nSqrtResult:%f\nTime:%d\n\n", num, result, iter_t);

	//割线法求根
	result = SecantSqrt(num, &iter_t);
	printf("割线法求根\nNum:%f\nSqrtResult:%f\nTime:%d\n\n", num, result, iter_t);
	
	system("Pause");
	return 0;
}

float NewtonSqrt(float c, int* iter_t) // c代表方根中的常数项
{
	// 迭代新值与旧值 0.5倍初值迭代
	float x_new = 0, x_pre = 0.5 * c;

	//统计迭代次数
	*iter_t = 1;

	//第一次迭代 迭代公式 x(n) = x(n-1) - f(x(n-1))/f'(x(n-1));
	x_new = x_pre - (x_pre * x_pre - c) / (2 * x_pre);

	//迭代到收敛为止 方根有唯一解 
	while (fabs(x_new - x_pre) >= err_limit && *iter_t <= iter_limit)
	{
		(*iter_t)++;
		x_pre = x_new;
		x_new = x_pre - (x_pre * x_pre - c) / (2 * x_pre);
	}

	//返回结果
	return x_new;
}

float DivisionSqrt(float c, int* iter_t)
{
	float x_left	= 0;
	float x_right	= c;
	float x_mid		= 0;

	//统计迭代次数
	*iter_t = 0;

	//迭代到收敛为止
	while (fabs(x_left - x_right) >= err_limit && *iter_t <= iter_limit)
	{
		(*iter_t)++;

		//折半查找
		x_mid = (x_left + x_right) / 2;
		
		if ((x_mid * x_mid - c) < 0) x_left  = x_mid;
		else                         x_right = x_mid;
	}

	//返回最新结果
	return x_mid;
}

float SecantSqrt(float c, int* iter_t)
{
	float x_new = c;
	float x_pre = 0;
	float x_cal = c;

	*iter_t = 0;

	//迭代到收敛为止
	while (fabs(x_new - x_pre) >= err_limit && *iter_t <= iter_limit)
	{
		(*iter_t)++;

		//迭代式
		//xn+1 = xn - ((xn-xn-1)/(f(xn)-f(xn-1)))*f(xn)
		x_cal = x_new - ((x_new - x_pre) / (x_new * x_new - x_pre * x_pre))*(x_new * x_new - c);
		x_pre = x_new;
		x_new = x_cal;
	}

	//返回最新结果
	return x_new;
}