newton.cpp

Newton法解方程

#include<iostream.h>
#include<math.h>
#include<iomanip.h>

const double e = 2.718281828;
void main()
{
	double p, p0, f0, f1, f2;
	int i = 0;
	p0 = 0.5;
	while(1){
		//i++;
		//p = p0 - (p0 * p0 - 10 * cos(p0)) / (2 * p0 + 10 * sin(p0));		
		//p = p0 - (pow(p0+1, 360) - 500*p0 - 1)/(360*pow(p0+1, 359) - 500);
		//p = p0 - ((p0*p0 - 2*p0*pow(e, (0-p0)) + pow(e, -2*p0))/(2*p0 - 2*pow(e, (0-p0)) + 2*p0*pow(e, (0-p0)) - 2*pow(e, -2*p0)));
		//p = p0 - (cos(p0+sqrt(2)) + p0*p0/2 + sqrt(2)*p0)/( (0-sin(p0+sqrt(2))) + p0 + sqrt(2));
		//p = p0 - (p0*p0*p0 - 3*p0*p0*pow(0.5, p0) + 3*p0*pow(0.25, p0) - pow(0.125, p0))/(3*p0*p0 - 6*p0*pow(0.5,p0) + 3*p0*p0*log(2)*pow(0.5,p0) + 3*pow(0.25, p0) - 3*p0*log(4)*pow(0.25, p0) + pow(0.125, p0)*log(8));
		//p = p0 - (pow(e, 6*p0) + 3*log(2)*log(2)*pow(e, 2*p0) - log(8)*pow(e, 4*p0) - log(2)*log(2)*log(2))/(6*pow(e,6*p0) + 6*log(2)*log(2)*pow(e, 2*p0) - 4*log(8)*pow(e, 4*p0));
		//f0 = p0*p0 - 2*p0*pow(e, 0-p0) + pow(e, -2*p0);
		//f1 = 2*p0 - 2*pow(e, (0-p0)) + 2*p0*pow(e, (0-p0)) - 2*pow(e, -2*p0);
		//f2 = 2 + 4*pow(e, 0-p0) - 2*p0*pow(e, 0-p0) + 4*pow(e, -2*p0);

		//f0 = cos(p0+sqrt(2)) + p0*p0/2 + sqrt(2)*p0;
		//f1 = (0-sin(p0+sqrt(2))) + p0 + sqrt(2);
		//f2 = 0 - cos(p0 + sqrt(2)) + 1;

		//f0 = p0*p0*p0 - 3*p0*p0*pow(0.5, p0) + 3*p0*pow(0.25, p0) - pow(0.125, p0);
		//f1 = 3*p0*p0 - 6*p0*pow(0.5,p0) + 3*p0*p0*log(2)*pow(0.5,p0) + 3*pow(0.25, p0) - 3*p0*log(4)*pow(0.25, p0) + pow(0.125, p0)*log(8);
		//f2 = 6*p0 - 6*pow(0.5, p0) + 12*p0*log(2)*pow(0.5, p0) - 3*p0*p0*log(2)*log(2)*pow(0.5, p0) - 6*log(4)*pow(0.25, p0) + 3*p0*log(4)*log(4)*pow(0.25, p0) - log(8)*log(8)*pow(0.125, p0);
		
		//f0 = pow(e, 6*p0) + 3*log(2)*log(2)*pow(e, 2*p0) - log(8)*pow(e, 4*p0) - log(2)*log(2)*log(2);
		//f1 = 6*pow(e,6*p0) + 6*log(2)*log(2)*pow(e, 2*p0) - 4*log(8)*pow(e, 4*p0);
		//f2 = 36*pow(e,6*p0) + 12*log(2)*log(2)*pow(e, 2*p0) - 16*log(8)*pow(e, 4*p0);
		
		//f0 = pow(e, 6*p0) + 1.441*pow(e, 2*p0) - 2.079*pow(e, 4*p0) - 0.333;
		//f1 = 6*pow(e, 6*p0) + 2.882*pow(e, 2*p0) - 8.316*pow(e, 4*p0);
		//f2 = 36*pow(e, 6*p0) + 5.764*pow(e, 2*p0) - 33.264*pow(e, 4*p0);
		
		f0 = 600*p0*p0*p0*p0 - 550*p0*p0*p0 + 200*p0*p0 - 20*p0 - 1;
		f1 = 2400*p0*p0*p0 - 1650*p0*p0 + 400*p0 - 20;
		//p = p0 - f0*f1/(f1*f1 - f0*f2);
		p = p0 - f0/f1;
		cout << i++ << " " << setprecision(7) << p0 << endl;
		if(fabs(p - p0) < pow(10, -4))
			break;
		p0 = p;
	}
	cout << p << endl;
}