「非線型方程式(後編)」(2009年4月号)のPDFファイルと本サポートページが3月号 のものになってましたので正しいものに差し替えました。
場所> | 正 | 誤 |
---|---|---|
p. 77右 (2)式分母 | nan+(n-1)b/2 | nan+(n-1)/2 |
テキストのPDFファイルをご覧下さい。
2009年3月号例4の出力に用いたもの。
/* Newton method */ #include <stdio.h> #include <math.h> #include <complex.h> #define IMAX 10 #define EPS 1.0e-12 // If |f(z_i)|< EPS then stop. double complex f(double complex z); double complex df(double complex z); int main(void) { int i; double absfz; double complex z,u,fz,dfz; z = -1.0 + 1.0 * I; i = 0; absfz = cabs(f(z)); printf("%2s %20s %20s %10s s\n","i","Re(z_i) ","Im(z_i) ","|f(z_i)|"); printf("%2d %20.15f %20.15f %10.3e\n",i,creal(z),cimag(z),absfz); while(absfz > EPS && i < IMAX) { fz = f(z); dfz = df(z); u = fz/dfz; z = z - u; i++; absfz = cabs(f(z)); printf("%2d %20.15f %20.15f %10.3e\n",i,creal(z),cimag(z),absfz); } return 0; } double complex f(double complex z) { double complex w; w=-5.0+z*(-2.0+z*z); return w; } double complex df(double complex z) { double complex w; w=-2.0 + 3.0 * z * z; return w; }
2009年4月号表4の出力に用いたもの。
// Weierstrass method with Aberth's radius #include <stdio.h> #include <math.h> #include <complex.h> #define PI 3.1415926535897932 #define NUMAX 30 // If nu >= NUMAX then stop. #define EPS 1.0e-09 //1.0e-12 // If max{|f(z_1)|,...,|f(z_n)|}< EPS then stop. #define N 5 // degree of the polynomial double complex f(double complex z, double complex *c); double newton(double x, double *b); double complex horner(double complex *a); double ff(double x, double *b); double dff(double x, double *b); int main(void) { int i,j,nu; double R,del; /* del:=|f(z_1)|+...+|f(z_n)| */ double complex z[N],z0[N],w; double complex c[N+1]={1.0,-3.0,9.0,-37.0,80.0,-50.0}; /* f(z)=c[0]*z^N+c[1]*z^{N-1}+...+c[N] */ double complex a[N+1]; double b[N+1]; for(i = 0;i <= N;i++) a[i] = c[i]; horner(a); /* Horner method */ b[0] = cabs(a[0]); for(i=1;i<=N;i++) b[i] = -cabs(a[i]); R = newton(100.0,b); /* 100.0 is the initial value */ printf(" R=%20.15f\n\n",R); /* for Table 4 */ R = 3.875; printf(" R=%20.15f\n\n",R); /* end Table 4 */ for(j = 0;j < N;j++){ z[j] = -c[1]/(c[0]*N)+R*cexp(I*PI*(2*j+0.5)/N); } del = 0.0; for(j = 0;j < N;j++){ if(del < cabs(f(z[j],c))) del = cabs(f(z[j],c)); } nu = 0; printf("%2d %10.3e\n",nu,del); for(j=0;jEPS && nu < NUMAX){ nu++; for(j = 0;j < N;j++){ z0[j] = z[j]; } for(i = 0;i < N;i++){ w = 1.0; for(j = 0;j < N;j++){ if (j == i) continue; w *= (z0[i] - z0[j]); } z[i] = z0[i] - f(z0[i],c)/(c[0]*w); } del = 0.0; for(j = 0;j < N;j++){ if(cabs(f(z[j],c)) > del) del = cabs(f(z[j],c)); } printf("%2d %10.3e\n",nu,del); for(j = 0;j < N;j++){ printf(" %d %20.15f %20.15f\n",j,creal(z[j]),cimag(z[j])); } } return 0; } double complex f(double complex z, double complex *c) // f(z) { int j; double complex w = c[0]; for(j = 1;j <= N;j++) w = c[j] + w*z; return w; } double ff(double x, double *b) { int j; double w = b[0]; for(j = 1;j <= N;j++) w = b[j]+w*x; return w; } double dff(double x, double *b) { int j; double w = N*b[0]; for(j = 1;j <= N-1;j++) w = (N-j)*b[j] + w*x; return w; } double complex horner(double complex *a) { int i,j; complex double gamma; gamma = -a[1]/(N * a[0]); for(i = 0;i <= N;i++){ for(j = 1;j <= N-i;j++){ a[j] = a[j] + gamma * a[j-1]; } } return 0; } double newton(double x, double *b) { int nu; nu=0; while(fabs(ff(x,b)) > EPS && nu < NUMAX){ x = x - ff(x,b)/dff(x,b); nu++; } return x; }
お気づきのことがありましたら、下記のメールアドレスにご一報ください。
メールアドレス 長田直樹 <osada "at" cis.twcu.ac.jp> (スパム防止のため@を"at"で置き換えています。)