ccsu 1027一元三次方程式の解(ニュートン反復)

rをf(x)=0のルートとし、x 0をr初期近似値として選択し、過点(x 0,f(x 0))を曲線y=f(x)の接線Lとし、Lの方程式をy=f(x 0)+f'(x 0)(x-x 0)とし、Lとx軸の交点の横座標x 1=x 0-f(x 0)/f'(x 0)を求め、x 1をrの一次近似値と呼ぶ.過点(x 1,f(x 1))は曲線y=f(x)の接線とし、この接線とx軸の交点の横座標x 2=x 1-f(x 1)/f'(x 1)を求め、x 2をrの二次近似値と呼ぶ.以上の手順を繰り返して、rの近似値シーケンスを得、x(n+1)=x(n)−f(x(n)/f′(x(n))をrのn+1次近似値と呼び、上式をニュートン反復式と呼ぶ.
非線形方程式f(x)=0を解くニュートン法は非線形方程式を線形化する近似法である.f(x)をx 0点付近でテイラー級数f(x)=f(x 0)+(x-x 0)f'(x 0)+(x-x 0)^2*f'(x 0)/2!+...その線形部分を取って、非線形方程式f(x)=0の近似方程式、すなわちテイラー展開の最初の2項として、f(x 0)+f′(x 0)(x-x 0)=0をf′(x 0)≠0とするとx 1=x 0-f(x 0)/f′(x 0)と解くと、ニュートン法の反復シーケンス:x(n+1)=x(n)−f(x(n))/f′(x(n))が得られる.

#include<iostream>
using namespace std;
double fabs(double x)
{
	if(x>0) return x;
	else return -x;
}
int main()
{
	double a,b,c,d;
	while (scanf("%lf%lf%lf%lf",&a,&b,&c,&d)!=EOF)
	{
		if(a==0&&b==0&&c==0&&d==0)  return 0;
		double x0=1,x,f;
		x=x0-(x0*x0*x0*a+x0*x0*b+c*x0+d)/(3*x0*x0*a+2*b*x0+c);
		f=x*x*x*a+x*x*b+c*x+d;
		while (fabs(f)-0>1e-8)
		{
			x0=x;
			x=x0-(x0*x0*x0*a+x0*x0*b+c*x0+d)/(3*x0*x0*a+2*b*x0+c);
	    	f=x*x*x*a+x*x*b+c*x+d;
		}
		printf("%.3lf
",x);
	}
	return 0;
}