http://www.
方程组
AX = b
如果A市方阵,当det(A)接近于0的时候,A的条件数很大,高斯消去法的一系列方法求解结果误差很大,甚至面目全非,这个时候使用矩阵平衡也几乎不起作用,只好用奇异值分解的方法
以下是从Numerical Recipes in C复制的代码,也许有人用的着
#define NR_END 1
#define FREE_ARG char*
void nrerror ( char error_text[] )
/* Numerical Recipes standard error handler */
{
fprintf ( stderr, "Numerical Recipes run-time error...\n" );
fprintf ( stderr, "%s\n", error_text );
fprintf ( stderr, "...now exiting to system...\n" );
exit ( 1 );
}
long double *fvector ( long nl, long nh )
/* allocate a long double vector with subscript range v[nl..nh] */
{
long double *v;
v = ( long double * ) malloc ( ( size_t ) ( ( nh - nl + 1 + NR_END ) * sizeof ( long double ) ) );
if ( !v ) nrerror ( "allocation failure in vector()" );
return v -nl + NR_END;
}
void free_fvector ( long double *v, long nl, long nh )
/* free a long double vector allocated with vector() */
{
free ( ( FREE_ARG ) ( v + nl - NR_END ) );
}
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
#define FMAX(a,b) ((a)>(b)? (a) : (b))
#define SQR(x) ((x)*(x))
#define IMIN(a,b) ((a) < (b) ?(a) : (b))
//---------------------------------------------------------------------------
void svdcmp ( long double **a, int m, int n, long double w[], long double **v )
{
long double pythag ( long double a, long double b );
int flag, i, its, j, jj, k, l, nm;
long double anorm, c, f, g, h, s, scale, x, y, z, *rv1;
rv1 = fvector ( 1, n );
.............................