cauchy.c

柯西方法解方程

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "math_util.h"
#define MAX_ITT_CM 5000

void cauchy( fun, dfun, x, n, tol )
long double (*fun)( long double * );
void (*dfun)( long double *, long double * );
long double *x, tol;
int n;
{
   long double dx, ndx, *grad, *h, *x0, *x1, f0, f1, alpha, err1;
   int i, j;
   long double line_search( long double(*)(), long double, long double *,
                             long double *, int, long double );
   
   if( tol <= 0.0 ) tol = 1.0e-4;
   if( n <= 0 ) return;

   grad = a1d_allo_dbl( n );
   h = a1d_allo_dbl( n );
   x0 = a1d_allo_dbl( n );
   
   vec_copy( x0, x, n );
   
   f0 = (*fun)( x0 );

   for( i = 0; i < MAX_ITT_CM; i++ ) {
      printf(" itt = %d, f = %Le\n",i,f0);

      (*dfun)( grad, x0 );
      err1 = enorm( grad, n );
      if( err1 <= tol ) {
         vec_copy( x, x0, n );
         a1d_free_dbl( x0 );
         a1d_free_dbl( h );
         a1d_free_dbl( grad );
         printf("***  Norm(grad f) <= eps\n");
         return;
      }
      for( j = 0; j < n; j++ )
         h[j] = -grad[j];
         
      alpha = line_search( fun, f0, x0, h, n, tol );
      
      ndx = 0.0;
      
      for( j = 0; j < n; j++ ) {
         dx = alpha*h[j];
         x0[j] += dx;
         ndx += dx*dx;
      }
      ndx = sqrt( ndx );
         
      f1 = (*fun)( x0 );
      
/*
      if( ndx <= tol/10 ) {
         vec_copy( x, x0, n );
         a1d_free_dbl( x0 );
         a1d_free_dbl( h );
         a1d_free_dbl( grad );
         printf("***  |x1-x0| <= 0.1*eps\n");
         return;
      }
*/   
      f0 = f1;   
   }
   printf("Error: Too many Cauchy iterations\n");
   vec_copy( x, x0, n );
   a1d_free_dbl( x0 );
   a1d_free_dbl( h );
   a1d_free_dbl( grad );
   return;   
}