muller.cpp

muller法解方程

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

double f(double x)
{
	return (600*x*x*x*x - 550*x*x*x + 200*x*x - 20*x - 1);
}

void main()
{
	double x2, x1, x0, h, h1, h2, d1, d2, d, p;
	int i = 3;
	
	x0 = 0.1;
	x1 = 0.2;
	x2 = 0.3;

	h1 = x1 - x0;
	h2 = x2 - x1;
	d1 = (f(x1) - f(x0))/h1;
	d2 = (f(x2) - f(x1))/h2;
	d = (d2 - d1)/(h2 + h1);

	double b, D, E;
	while(1)
	{
		b = d2 + h2*d;
		D = sqrt(b*b - 4*f(x2)*d);
		if(fabs(b-D) < fabs(b+D)) E = b+D;
		else E = b - D;
		h = -2*f(x2)/E;
		p = x2 + h;
		if(fabs(h) < 0.0001){
			cout << p << endl;
			break;
		}
		else{
			cout << p << endl;
			x0 = x1;
			x1 = x2;
			x2 = p;
			h1 = x1 - x0;
			h2 = x2 - x1;
			d1 = (f(x1) - f(x0))/h1;
			d2 = (f(x2) - f(x1))/h2;
			d = (d2 - d1) / (h2 + h1);
			i ++;
		}
	}
}