「非線型方程式(後編)」(2009年4月号)

「非線型方程式(後編)」(2009年4月号)のPDFファイルと本サポートページが３月号 のものになってましたので正しいものに差し替えました。

正誤表

テキストのPDFファイルをご覧下さい。

プログラム

複素Newton法のプログラム

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;
}

Weierstrass法のプログラム

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;j EPS && 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"で置き換えています。)
Last modified: Wed Sep 23 10:30:47 JST 2009