atl.rar_图形图象

/* * MULLER'S ALGORITHM 2.8 * * To find a solution to f(x) = 0 given three approximations x0, x1 * and x2: * * INPUT: x0,x1,x2; tolerance TOL; maximum number of iterations NO. * * OUTPUT: approximate solution p or message of failure. * * This implementation allows for a switch to complex arithmetic. * The coefficients are stored in the vector A, so the dimension * of A may have to be changed. */ #include<stdio.h> #include<math.h> #define ZERO 1.0E-20 #define pi 3.141592654 #define true 1 #define false 0 double absval(double); void CADD(double, double, double, double, double *, double *); void CMULT(double, double, double, double, double *, double *); void CDIV(double, double, double, double, double *, double *); double CABS(double, double); void CSQRT(double, double, double *, double *); void CSUB(double, double, double, double, double *, double *); void INPUT(int *, double *, double *, double *, int *, int *); void OUTPUT(FILE **); main() { double A[50],ZR[4],ZI[4],GR[4],GI[4],F[4],X[4],CH1R[3],CH1I[3],H[3]; double CDEL1R[2],CDEL1I[2],DEL1[2]; double CDELR,CDELI,CBR,CBI,CDR,CDI,CER,CEI,DEL,B,D,E,TOL,QR,QI,ER,EI,FR,FI; int N,M,I,J,K,ISW,KK,OK; FILE *OUP[1]; printf("This is Mullers Method.\n"); INPUT(&OK, A, X, &TOL, &M, &N); if (OK) { OUTPUT(OUP); /* Evaluate F using Horner's Method and save in the vector F */ for (I=1; I<=3; I++) { F[I-1] = A[N-1]; for (J=2; J<=N; J++) { K = N-J+1; F[I-1] = F[I-1]*X[I-1]+A[K-1]; } } /* Variable ISW is used to note a switch to complex arithmetic ISW=0 means real arithmetic, and ISW=1 means complex arithmetic */ ISW = 0; /* STEP 1 */ H[0] = X[1]-X[0]; H[1] = X[2]-X[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]); I = 3; OK = true; /* STEP 2 */ while ((I <= M) && OK) { /* STEPS 3-7 for real arithmetic */ if (ISW == 0) { /* STEP 3 */ B = DEL1[1]+H[1]*DEL; D = B*B-4.0*F[2]*DEL; /* test to see if need complex arithmetic */ if (D >= 0.0) { /* real arithmetic - test to see if straight line */ if (absval(DEL) <= ZERO) { /* straight line - test if horizontal line */ if (absval(DEL1[1]) <= ZERO) { printf("Horizontal Line\n"); OK = false; } else { /* straight line but not horizontal */ X[3] = (F[2]-DEL1[1]*X[2])/DEL1[1]; H[2] = X[3]-X[2]; } } else { /* not straight line */ D = sqrt(D); /* STEP 4 */ E = B+D; if (absval(B-D) > absval(E)) E = B-D; /* STEP 5 */ H[2] = -2.0*F[2]/E; X[3] = X[2]+H[2]; } if (OK) { /* evaluate f(x(I))=F[3] by Horner's method */ F[3] = A[N-1]; for (J=2; J<=N; J++) { K = N-J+1; F[3] = F[3]*X[3]+A[K-1]; } fprintf(*OUP, "%d %f %f\n",I,X[3],F[3]); /* STEP 6 */ if (absval(H[2]) < TOL) { fprintf(*OUP, "\nMethod Succeeds\n"); fprintf(*OUP, "Approximation is within %.10e\n",TOL); fprintf(*OUP, "in %d iterations\n", I); OK = false; } else { /* STEP 7 */ 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 { /* switch to complex arithmetic */ ISW = 1; for (J=1; J<=3; J++) { ZR[J-1] = X[J-1]; ZI[J-1] = 0.0; GR[J-1] = F[J-1]; GI[J-1] = 0.0; } for (J=1; J<=2; J++) { CDEL1R[J-1] = DEL1[J-1]; CDEL1I[J-1] = 0.0; CH1R[J-1] = H[J-1]; CH1I[J-1] = 0.0; } CDELR = DEL; CDELI = 0.0; } } if ((ISW == 1) && OK) { /* STEPS 3-7 complex arithmetic */ /* test if straight line */ if (CABS(CDELR,CDELI) <= ZERO) { /* straight line - test if horizontal line */ if (CABS(CDEL1R[0],CDEL1I[0]) <= ZERO) { printf("horizontal line - complex\n"); OK = false; } else { /* straight line but not horizontal */ printf("line - not horizontal\n"); CMULT(CDEL1R[1],CDEL1I[1],ZR[2],ZI[2],&ER,&EI); CSUB(GR[2],GI[2],ER,EI,&FR,&FI); CDIV(FR,FI,CDEL1R[1],CDEL1I[1],&ZR[3],&ZI[3]); CSUB(ZR[3],ZI[3],ZR[2],ZI[2],&CH1R[2],&CH1I[2]); } } else { /* not straight line */ /* STEP 3 */ CMULT(CH1R[1],CH1I[1],CDELR,CDELI,&ER,&EI); CADD(CDEL1R[1],CDEL1I[1],ER,EI,&CBR,&CBI); CMULT(GR[2],GI[2],CDELR,CDELI,&ER,&EI); CMULT(ER,EI,4.0,0.0,&FR,&FI); QR = CBR; QI = CBI; CMULT(CBR,CBI,QR,QI,&ER,&EI); CSUB(ER,EI,FR,FI,&CDR,&CDI); CSQRT(CDR,CDI,&FR,&FI); /* STEP 4 */ CDR = FR; CDI = FI; CADD(CBR,CBI,CDR,CDI,&CER,&CEI); CSUB(CBR,CBI,CDR,CDI,&FR,&FI); if (CABS(FR,FI) > CABS(CER,CEI)) CSUB(CBR,CBI,CDR,CDI,&CER,&CEI); /* STEP 5 */ CDIV(GR[2],GI[2],CER,CEI,&ER,&EI); CMULT(ER,EI,-2.0,0.0,&CH1R[2],&CH1I[2]); CADD(ZR[2],ZI[2],CH1R[2],CH1I[2],&ZR[3],&ZI[3]); } if (OK) { /* evaluate f(x(i))=f(3) by Horner's method */ GR[3] = A[N-1]; GI[3] = 0.0; for (J=2; J<=N; J++) { K = N-J+1; CMULT(GR[3],GI[3],ZR[3],ZI[3],&ER,&EI); GR[3] = ER+A[K-1]; GI[3] = EI; } fprintf(*OUP, "%4d %14.8f %14.8f %14.8f %14.8f\n", I, ZR[3], ZI[3], GR[3], GI[3]); /* STEP 6*/ if (CABS(CH1R[2],CH1I[2]) <= TOL) { fprintf(*OUP, "\nMethod Succeeds\n"); fprintf(*OUP, "Approximation is within %.8e\n", TOL); fprintf(*OUP, "in %d iterations\n", I); OK = false; } else { /* STEP 7 */ for (J=1; J<=2; J++) { CH1R[J-1] = CH1R[J]; CH1I[J-1] = CH1I[J]; ZR[J-1] = ZR[J]; ZI[J-1] = ZI[J]; GR[J-1] = GR[J]; GI[J-1] = GI[J]; } ZR[2] = ZR[3]; ZI[2] = ZI[3]; GR[2] = GR[3]; GI[2] = GI[3]; CDEL1R[0] = CDEL1R[1