牛顿下山法求解非线性方程（组）（C实现）

1、算法描述

（1）符号说明与基本假设

对于非线性方程组：

                                                       （1）

引入向量：

可将（1）式改写为

                                                       （2）

通常考虑方程（2）只有唯一解的情形。

（2）牛顿下山算法

引入下山因子λ，产生一下计算格式：

下山因子λ一般依次取1、1/2、1/4、1/8、……

其中

计算步骤为：

2、C语言实现

newton.h头文件：

#ifndef __NEWTON_H__
#define __NEWTON_H__
// 牛顿下山法求解非线性方程（组）

int newton( double** X0
	, int n
	, double lmada
	, double eps_x
	, double eps_f
	, void   (*f)( double** X, int n)
	, void   (*df)(double** X, int n)
	);

#endif

newton.c文件

#include "newton.h"
#include <stdlib.h>
#include <string.h>

// 计算矩阵的逆
static int inverse( double** dfX, int n )
{
	int i, j,k, p;
	double maxV, tmp;
	double* A = *dfX;
	double* T = (double*)malloc(sizeof(double) * n * n);
	double* tArr = (double*)malloc(sizeof(double) * n);
	int row_size = sizeof(double) * n;

	memset( T, 0, sizeof(double)*n*n);
	for ( i = 0; i < n; ++i )
	{
		T[i*n+i] = 1.0;
	}

	for ( j = 0; j < n; ++j )
	{
		p = j;
		tmp = A[j*n+j];
		maxV = (tmp>=0)?(tmp):(-tmp);
		for ( i = j +1; i < n; ++i )
		{
			tmp = A[i*n+j];
			if ( tmp < 0 ) tmp = -tmp;
			if ( maxV < tmp )
			{
				p = i;
				maxV = tmp;
			}
		}

		if ( maxV < 1e-20 )
		{
			return -1;
		}

		if ( j != p )
		{
			memcpy(  tArr, A+j*n, row_size);
			memcpy( A+j*n, A+p*n, row_size);
			memcpy( A+p*n,  tArr, row_size);
			memcpy(  tArr, T+j*n, row_size);
			memcpy( T+j*n, T+p*n, row_size);
			memcpy( T+p*n,  tArr, row_size);
		}

		tmp = A[j*n+j];
		for ( i = j; i < n; ++i ) A[j*n+i] /= tmp;
		for ( i = 0; i < n; ++i ) T[j*n+i] /= tmp;
		for ( i = 0; i < n; ++i )
		{
			if ( i != j )
			{
				tmp = A[i*n+j];
				for ( k = j; k < n; ++k )
					A[i*n+k] -= tmp * A[j*n+k];
				for ( k = 0; k < n; ++k )
					T[i*n+k] -= tmp * T[j*n+k];
			}
		}
	}
	memcpy( A, T, row_size * n );
	free( T );
	free( tArr );

	return 0;
}

// 计算步长dx
static void calc_dx( double** dx
	               , double*  df
				   , double*  dfx
				   , double   lamda
				   , int      n
				   )
{
	int i, j;
	double* x = *dx;
	memset( x, 0, sizeof(double) * n);
	for ( i = 0; i < n; ++i )
	{
		for ( j = 0; j < n; ++j )
		{
			x[i] += -lamda * df[j] * dfx[i*n+j];
		}
	}
}

// 计算向量的无穷范数
static double norm_inf( double* A, int n )
{
	int i;
	double t = A[0];
	double ret = t;

	if ( t < 0 )
	{
		ret = -t;
	}

	for ( i = 1; i < n; ++i )
	{
		t = A[i];
		if ( t < 0 ) t = -t;
		if ( ret < t ) ret = t;
	}

	return ret;
}

// 牛顿下山法求解非线性方程组，求解成功返回0，失败返回-1
int  newton  ( double** X0                     // 迭代起始点
	          , int n                           // 方程组维数
			  , double lamda                    // 起始下山因子
			  , double eps_x                    // 阈值
			  , double eps_f                    // 阈值
			  , void   (*f)( double** X, int n) // 带求解非线性方程组函数
			  , void   (*df)(double** X, int n) // 带求解非线性方程组的导函数（Jacobi矩阵）
			  )
{
	int i, ret = 0;
	int row_size = sizeof(double) * n;
	double*   X = *X0;
	double*  dx = (double*)malloc( row_size );
	double*  fX = (double*)malloc( row_size );
	double* dfX = (double*)malloc( n * row_size );
	double max_f, max_f1;

	memcpy( fX, X, row_size );
	f ( &fX, n );
	for ( ;; )
	{
		memcpy( dfX, X, row_size);
		df( &dfX, n );

		ret = inverse( &dfX, n ); // 计算逆
		
		if ( ret < 0 ) // Jacobi矩阵不可逆
		{
			goto end;
		}
		calc_dx( &dx, fX, dfX, lamda, n); // 计算步长
		
		max_f = norm_inf( fX, n );
		for ( i = 0; i < n; ++i )
		{
			X[i] += dx[i];
		}
		
		memcpy( fX, X, row_size);
		f( &fX, n );
		max_f1 = norm_inf( fX, n );
		
		if ( max_f1 < max_f )
		{
			if ( norm_inf( dx, n ) < eps_x )
			{
				break;
			}
			else
			{
				continue;
			}
		}
		else
		{
			if ( max_f1 < eps_f )
			{
				break;
			}
			else
			{
				lamda /= 2.0;
			}
		}
	}

end:
	free( dx );
	free( fX );
	free( dfX);

	return ret;
}

3、测试

考虑用牛顿下山法求解以下非线性方程组：

求解以上方程的主程序：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#include "newton.h"
#define  MATH_E   2.7182818285
#define  MATH_PI  3.1415926536


void f( double** X, int n )
{
	assert( n == 3 );
	double* A = *X;
	double x = A[0];
	double y = A[1];
	double z = A[2];

	A[0] = 3*x-cos(y*z)-0.5;
	A[1] = x*x-81*(y+0.1)*(y+0.1)+sin(z)+1.06;
	A[2] = pow( MATH_E, -x*y)+20*z + 10*MATH_PI/3.0-1;
}

void df( double** X, int n )
{
	assert( n == 3 );
	double* A = *X;
	double x = A[0];
	double y = A[1];
	double z = A[2];

	A[0] = 3.0;
	A[1] = z * sin(y*z);
	A[2] = y * sin(y*z);
	A[3] = 2*x;
	A[4] = -162.0*(y+0.1);
	A[5] = cos(z);
	A[6] = -y * pow( MATH_E, -x*y);
	A[7] = -x * pow( MATH_E, -x*y);
	A[8] = 20.0;
}

int main()
{
	int n = 3;
	double* X = (double*)malloc(sizeof(double)*n);
	X[0] = 1.0;
	X[1] = 1.0;
	X[2] = 1.0;

	double eps_x = 1e-14;
	double eps_f = eps_x;
	double lamda = 1.0;

	newton( &X, n, lamda, eps_x, eps_f, f, df);
	printf("%f\t%f\t%f", X[0], X[1], X[2]);

	free( X );
	return 1;
}

计算结果：

x = 0.500000

y = -0.000000

z = -0.523599