牛顿下山法求解非线性方程(组)(C实现)
2012-04-17 18:25
267 查看
1、算法描述
(1)符号说明与基本假设
对于非线性方程组:
(1)
引入向量:
可将(1)式改写为
(2)
通常考虑方程(2)只有唯一解的情形。
(2)牛顿下山算法
引入下山因子λ,产生一下计算格式:
下山因子λ一般依次取1、1/2、1/4、1/8、……
其中
计算步骤为:
2、C语言实现
newton.h头文件:
newton.c文件
3、测试
考虑用牛顿下山法求解以下非线性方程组:
求解以上方程的主程序:
计算结果:
x = 0.500000
y = -0.000000
z = -0.523599
(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
相关文章推荐
- 非线性方程组求解算法的C++实现
- matlab实现牛顿迭代法求解非线性方程组
- matlab实现牛顿迭代法求解非线性方程组
- 使用C语言实现最小生成树求解的简单方法
- 类动态规划求解较小规模的最大团问题(Python实现)
- 怎么实现每隔3个数字换行,求解?
- C语言实现Hanoi(汉诺塔)问题的递归求解
- 蚁群算法java实现以及TSP问题蚁群算法求解
- 同 用Java实现天平称球问题的自动求解
- 高斯消元法求解线性方程(c实现)
- 求解数独的C++实现
- 递归求解数组的和 Java实现
- 用栈实现后缀表达式求解问题
- Python求解非线性方程
- 递归实现最大值求解
- C++实现二分法求解方程
- 关于递归的总结——汉诺塔、素因数的求解(Python实现)
- Python实现Json结构对比的小工具兼谈编程求解问题
- 算术表达式求解java实现
- X的N次方求解——pow(x,n)实现