// Muller의 방법
#include <complex.h> // Borland-C에서는 복소수 클래스가 지원됨.
#include <stdio.h>
#include <math.h>
#include <iostream.h>
#define ZERO 1.0E-20
#define pi 3.141592654
#define true 1
#define false 0
void main()
{
complex Z_[4], G_[4], CH1_[3], CDEL1_[2], CDEL_, CB_, CD_, CE_, Q_, E_, F_;
// 복소수 변수는 변수명 뒤에 _ 를 붙여 표시
double A[50], F[4], X[4], H[3];
double DEL1[2];
double DEL, B, D, E;
int degree, i, j, KK, OK, isComplex;
double tolerance = 1.e-4; // 오차의 한계
int num = 40; // 계산의 횟수
FILE *display;
// 출력위치 지정
cout << "Select Output :" << endl;
cout << " 1 : screen (default)" << endl;
cout << " 2 : file (->muller.txt)" << endl;
cin >> i;
if(i == 2) display = fopen("muller.txt","w");
else display = stdout;
// 2번을 선택하면 파일로 저장, 나머지 경우는 화면 출력
fprintf(display, "Muller's Method.\n");
OK = false;
while (!OK) {
cout << "Input the degree of the polynomial" << endl;
cin >> degree; // 차수를 입력받는다.
if (degree++ > 0) {
OK = true;
cout << "\nInput the coefficients of the " \
"polynomial in ascending order\n" \
"of exponent at the prompt" << endl;
for (i=0; i<degree ; i++) {
cout << "Input A(" << i << ")\n";
cin >> A[i]; // 계수를 입력받는다.
}
} else cout << "n must be a positive integer" << endl;
}
if (A[degree-1] ==0) { // 최고차항 계수가 0일 경우
cout << "Leading coefficient is 0 - error in input" << endl;
OK = false;
}
if (degree == 2) { // 1차식인 경우
printf("Polynomial is linear: root is %11.8f", -A[0]/A[1]);
OK = false;
}
if (OK) {
OK = true;
//계산할 세 점을 받아들인다.
cout << "Input the first of three starting values" << endl;
scanf("%lf", &X[0]);
cout << "Input the second of three starting values" << endl;
scanf("%lf", &X[1]);
cout << "Input the third starting value" << endl;
scanf("%lf", &X[2]);
}
if (OK) {
// 주어진 세점을 이용해서 세 점을 지나는 함수의 값을 F[0]~F[2] 에 저장
for (i=0; i<3; i++) { // Horner's Method 를 사용
F[i] = A[degree-1];
for (j=1; j<degree; j++)
F[i] = F[i] * X[i] + A[degree-j-1];
}
// isComplex = false 이면 실수 계산이고, true 이면 복소수 계산이다.
isComplex = false;
H[0] = X[1] - X[0]; // H[0]은 처음 두점 사이의 거리
H[1] = X[2] - X[1]; // H[1]은 다음 두점 사이의 거리
DEL1[0] = (F[1] - F[0]) / H[0]; // 처음 두 점 사이의 기울기
DEL1[1] = (F[2] - F[1]) / H[1]; // 나중 두 점 사이의 기울기
DEL = (DEL1[1] - DEL1[0]) / (H[1] + H[0]);
// DEL은 세 점에서 구한 2차항의 계수(== a)
i = 3;
OK = true;
while ((i++ <= num) && OK) {
if (isComplex == false) { // 실수 계산일때
B = DEL1[1] + H[1] * DEL; // B는 세 점에서 구한 1차항의 계수(==b)
D = B * B - 4.0 * F[2] * DEL; // 판별식 (B==b, F[2]==c, DEL==a)
if (D >= 0.0) { // 판별식 - 실근일 경우
// 그래프가 1차 직선인 경우, 즉 최고차항의 계수가 0인 경우
if (fabs(DEL) <= ZERO) {
// 나중 두 점 사이의 기울기가 0일 경우
if (fabs(DEL1[1]) <= ZERO) {
printf("Horizontal Line\n");
OK = false;
} else { // 기울기가 0이 아닐경우
X[3] = (F[2]-DEL1[1]*X[2])/DEL1[1];
H[2] = X[3]-X[2];
}
}
else { // 그래프가 1차 직선이 아닌 경우
D = sqrt(D);
E = B+D;
if (fabs(B-D) > fabs(E)) E = B-D;
H[2] = -2.0*F[2]/E;
X[3] = X[2]+H[2];
}
if (OK) {
F[3] = A[degree-1]; // Horner's method를 사용해서 f(3)을 구한다.
for (j=1; j<degree; j++) F[3] = F[3]*X[3]+A[degree-j-1];
fprintf(display, "x(%d) -> %16.14f\n", i, X[3]);
fprintf(display, "F(x(%d)) -> %16.14f\n\n", i, F[3]);
// 허용오차보다 작을 경우 계산을 끝낸다.
if (fabs(H[2]) < tolerance) {
fprintf(display, "Results are real Number.\n");
fprintf(display, "\nMethod Succeeds");
fprintf(display, "\nApproximation -> %16.14e", X[3]);
fprintf(display, "\nTolerance -> %.10e", tolerance);
fprintf(display, "\nIterations -> %d\n", i);
OK = false;
}
else {
for (j=1; j<=2; j++) {
H[j-1] = H[j];
X[j-1] = X[j];
F[j-1] = F[j];
}
X[2] = X[3];
F[2] = F[3];
DEL1[0] = DEL1[1];
DEL1[1] = (F[2]-F[1])/H[1];
DEL = (DEL1[1]-DEL1[0])/(H[1]+H[0]);
}
}
}
else { // 판별식에서 허근이 나올 경우
isComplex = 1;
for (j=1; j<=3; j++) {
Z_[j-1] = complex(X[j-1], 0.0);
G_[j-1] = complex(F[j-1], 0.0);
}
for (j=1; j<=2; j++) {
CDEL1_[j-1] = complex(DEL1[j-1], 0.0);
CH1_[j-1] = complex(H[j-1], 0.0);
}
CDEL_ = complex(DEL, 0.0);
}
}
if ((isComplex == true) && OK) { // 허수 계산일 경우
// 최고차항의 계수가 0인 경우(그래프가 1차 직선)
if ( abs(CDEL_) <= ZERO) {
// 기울기가 0인 경우
if ( abs(CDEL1_[0]) <= ZERO) {
printf("horizontal line - complex\n");
OK = false;
}
else { // 기울기가 0이 아닌 경우
printf("line - not horizontal\n");
E_ = CDEL1_[1] * Z_[2];
F_ = G_[2] - E_;
Z_[3] = F_ / CDEL1_[1];
CH1_[2] = Z_[3] - Z_[2];
}
}
else { // 최고차항이 0이 아닌 경우(그래프가 2차 곡선)
E_ = CH1_[1] * CDEL_;
CB_ = CDEL1_[1] + E_;
E_ = G_[2] * CDEL_;
F_ = E_ * complex(4.0, 0.0);
Q_ = CB_;
E_ = CB_ * Q_;
CD_ = E_ - F_;
F_ = sqrt(CD_);
CD_ = F_;
CE_ = CB_ + CD_;
F_ = CB_ - CD_;
if (abs(F_) > abs(CE_)) CE_ = CB_ - CD_;
E_ = G_[2] / CE_;
CH1_[2] = E_ * complex(-2.0, 0.0);
Z_[3] = Z_[2] + CH1_[2];
}
if (OK) {
G_[3] = complex(A[degree-1],0.0);
// Horner's method를 사용해서 복소수 f(3)을 구함.
for (j=1; j<degree; j++) {
E_ = G_[3] * Z_[3];
G_[3] = complex(real(E_)+A[degree-j-1], imag(E_));
}
fprintf(display, "x(%d) -> %16.14f %+16.14f i\n", i, real(Z_[3]), imag(Z_[3]));
fprintf(display, "F(x(%d)) -> %16.14f %+16.14f i\n\n",i, real(G_[3]), imag(G_[3]));
// 허용오차보다 작을 경우 계산을 끝낸다.
if (abs(CH1_[2]) <= tolerance) {
fprintf(display, "Results are complex Number\n");
fprintf(display, "\nMethod Succeeds\n");
fprintf(display, "\nApproximation -> %16.14f %+16.14f i", real(Z_[3]), imag(Z_[3]));
fprintf(display, "\nTolerance -> %.10e", tolerance);
fprintf(display, "\nIterations -> %d\n", i);
OK = false;
}
else {
for (j=1; j<=2; j++) {
CH1_[j-1] = CH1_[j];
Z_[j-1] = Z_[j];
G_[j-1] = G_[j];
}
Z_[2] = Z_[3];
G_[2] = G_[3];
CDEL1_[0] = CDEL1_[1];
E_ = G_[2] - G_[1];
CDEL1_[1] = E_ / CH1_[1];
E_ = CDEL1_[1] - CDEL1_[0];
F_ = CH1_[1] + CH1_[0];
CDEL_ = E_ / F_;
}
}
}
}
// 계산횟수만큼 계산해도 허용오차보다 클 경우
if (OK) fprintf(display, "Method failed\n");
fclose(display);
}
}