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

#include

// 计算矩阵的逆

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

#include

#include

#include

#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