c语言求解病态方程组,病态方程组求解方法的C++实现（支持复系数方程组）

最新推荐文章于 2023-11-30 11:45:21 发布

weixin_39939918

最新推荐文章于 2023-11-30 11:45:21 发布

阅读量294

点赞数

文章标签： c语言求解病态方程组

头文件：

* This program is free software; you can redistribute it and/or modify it

* under the terms of the GNU General Public License as published by the

* Free Software Foundation, either version 2 or any later version.

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions are met:

* 1. Redistributions of source code must retain the above copyright notice,

* this list of conditions and the following disclaimer.

* 2. Redistributions in binary form must reproduce the above copyright

* notice, this list of conditions and the following disclaimer in the

* documentation and/or other materials provided with the distribution.

* This program is distributed in the hope that it will be useful, but WITHOUT

* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for

* more details. A copy of the GNU General Public License is available at:

* http://www.fsf.org/licensing/licenses

/*****************************************************************************

* linequs3.h

* Function template for solving Rank Defect linear equations.

* For a m-by-n coefficient matrix A and m-by-1 constant vector b, if A is a

* Rank Deficient Matrix, then the conventional methods will fail to solve

* Ax = b. Three methods for solving such problem are provided in this file,

* they are: Truncated SVD, Dampted SVD and Tikhonov Regularization.

* These methods adapt to both REAL or COMPLEX linear equations.

* Zhang Ming, 2010-07 (revised 2010-12), Xi'an Jiaotong University.

*****************************************************************************/

#ifndef RANKDEFLINEQUS_H

#define RANKDEFLINEQUS_H

#include

namespace splab

{

template

Vector tsvd( const Matrix&, const Vector&,

Real tol=Real(-1.0) );

template

Vector > tsvd( const Matrix >&,

const Vector >&, Type tol=Type(-1.0) );

template

Vector dsvd( const Matrix&, const Vector&, Real& );

template

Vector > dsvd( const Matrix >&,

const Vector >&, Type& );

template

Vector tikhonov( const Matrix&, const Vector&, Real& );

template

Vector > tikhonov( const Matrix >&,

const Vector >&, Type& );

#include

}

// namespace splab

#endif

// RANKDEFLINEQUS_H

实现文件：

* This program is free software; you can redistribute it and/or modify it

* under the terms of the GNU General Public License as published by the

* Free Software Foundation, either version 2 or any later version.

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions are met:

* 1. Redistributions of source code must retain the above copyright notice,

* this list of conditions and the following disclaimer.

* 2. Redistributions in binary form must reproduce the above copyright

* notice, this list of conditions and the following disclaimer in the

* documentation and/or other materials provided with the distribution.

* This program is distributed in the hope that it will be useful, but WITHOUT

* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for

* more details. A copy of the GNU General Public License is available at:

* http://www.fsf.org/licensing/licenses

/*****************************************************************************

* linequs3-impl.h

* Implementation for Rank Defect linear equations.

* Zhang Ming, 2010-07 (revised 2010-12), Xi'an Jiaotong University.

*****************************************************************************/

/**

* Rank defect linear equationequations solution by Truncated SVD.

* A ---> The m-by-n coefficient matrix(Rank Defect);

* b ---> The n-by-1 right-hand side vector;

* x ---> The n-by-1 solution vector.

template

Vector tsvd( const Matrix &A, const Vector &b, Real tol )

{

int m = A.rows(),

n = A.cols();

assert( m == b.size() );

SVD svd;

svd.dec( A );

Matrix U = svd.getU();

Matrix V = svd.getV();

Vector s = svd.getSV();

Vector x(n);

int r = 0;

if( tol <= 0 )

tol = max( m, n ) * s[0] * EPS;

for( int i=0; i

if( s[i] >= tol )

r++;

// y = U^T * b

Vector y(r);

for( int i=0; i

for( int j=0; j

y[i] += U[j][i]*b[j];

// y = y / s

for( int i=0; i

y[i] /= s[i];

// x = V * y

for( int i=0; i

for( int j=0; j

x[i] += V[i][j]*y[j];

// for( int i=0; i

// for( int j=0; j

// {

// Real sum = 0;

// for( int k=0; k

// sum += U[k][j]*b[k];

// x[i] += sum * V[i][j] / s[j];

// }

return x;

}

/**

* Rank defect complex linear equationequations solution by Truncated SVD.

* A ---> The m-by-n coefficient matrix(Rank Defect);

* b ---> The n-by-1 right-hand side vector;

* x ---> The n-by-1 solution vector.

template

Vector > tsvd( const Matrix > &A,

const Vector > &b,

Type tol )

{

int m = A.rows(),

n = A.cols();

assert( m == b.size() );

CSVD svd;

svd.dec( A );

Matrix > U = svd.getU();

Matrix > V = svd.getV();

Vector s = svd.getSV();

Vector > x(n);

int r = 0;

if( tol <= 0 )

tol = max( m, n ) * s[0] * EPS;

for( int i=0; i

if( s[i] >= tol )

r++;

// y = U^H * b

Vector > y(r);

for( int i=0; i

for( int j=0; j

y[i] += conj(U[j][i])*b[j];

// y = y / s

for( int i=0; i

y[i] /= s[i];

// x = V * y

for( int i=0; i

for( int j=0; j

x[i] += V[i][j]*y[j];

// for( int i=0; i

// for( int j=0; j

// {

// complex sum = 0;

// for( int k=0; k

// sum += conj(U[k][j])*b[k];

// x[i] += sum * V[i][j] / s[j];

// }

return x;

}

/**

* Rank defect linear equationequations solution by Dampted SVD.

* A ---> The m-by-n(m>n) coefficient matrix(Rank Defect);

* b ---> The n-by-1 right-hand side vector;

* sigma : dampted factor;

* x ---> The n-by-1 solution vector.

template

Vector dsvd( const Matrix &A, const Vector &b,

Real &sigma )

{

int m = A.rows(),

n = A.cols();

assert( m == b.size() );

SVD svd;

svd.dec( A );

Matrix U = svd.getU();

Matrix V = svd.getV();

Vector s = svd.getSV();

Vector x(n);

int p = s.size();

s += sigma;

// y = U^T * b

Vector y(p);

for( int i=0; i

for( int j=0; j

y[i] += U[j][i]*b[j];

// y = y / s

for( int i=0; i

y[i] /= s[i];

// x = V * y

for( int i=0; i

for( int j=0; j

x[i] += V[i][j]*y[j];

// for( int i=0; i

// for( int j=0; j

// {

// Real sum = 0;

// for( int k=0; k

// sum += U[k][j]*b[k];

// x[i] += sum * V[i][j] / s[j];

// }

return x;

}

/**

* Rank defect complex linear equationequations solution by Dampted SVD.

* A ---> The m-by-n(m>n) coefficient matrix(Rank Defect);

* b ---> The n-by-1 right-hand side vector;

* sigma : dampted factor;

* x ---> The n-by-1 solution vector.

template

Vector > dsvd( const Matrix > &A,

const Vector > &b,

Type &sigma )

{

int m = A.rows(),

n = A.cols();

assert( m == b.size() );

CSVD svd;

svd.dec( A );

Matrix > U = svd.getU();

Matrix > V = svd.getV();

Vector s = svd.getSV();

Vector > x(n);

int p = s.size();

s += sigma;

// y = U^H * b

Vector > y(p);

for( int i=0; i

for( int j=0; j

y[i] += conj(U[j][i])*b[j];

// y = y / s

for( int i=0; i

y[i] /= s[i];

// x = V * y

for( int i=0; i

for( int j=0; j

x[i] += V[i][j]*y[j];

// for( int i=0; i

// for( int j=0; j

// {

// complex sum = 0;

// for( int k=0; k

// sum += conj(U[k][j])*b[k];

// x[i] += sum * V[i][j] / s[j];

// }

return x;

}

/**

* Rank defect linear equationequations solution by Tikhonov Regularization.

* A ---> The m-by-n(m>n) coefficient matrix(Rank Defect);

* b ---> The n-by-1 right-hand side vector;

* alpha : regularization factor;

* x ---> The n-by-1 solution vector.

template

Vector tikhonov( const Matrix &A, const Vector &b,

Real &alpha )

{

int m = A.rows(),

n = A.cols();

assert( m == b.size() );

SVD svd;

svd.dec( A );

Matrix U = svd.getU();

Matrix V = svd.getV();

Vector s = svd.getSV();

Vector x(n);

int p = s.size();

Real alpha2 = alpha*alpha;

for( int i=0; i

s[i] += alpha2/s[i];

// y = U^T * b

Vector y(p);

for( int i=0; i

for( int j=0; j

y[i] += U[j][i]*b[j];

// y = y / s

for( int i=0; i

y[i] /= s[i];

// x = V * y

for( int i=0; i

for( int j=0; j

x[i] += V[i][j]*y[j];

// for( int i=0; i

// for( int j=0; j

// {

// Real sum = 0;

// for( int k=0; k

// sum += U[k][j]*b[k];

// x[i] += sum * V[i][j] / s[j];

// }

return x;

}

/**

* Rank defect complex linear equationequations solution by

* Tikhonov Regularization.

* A ---> The m-by-n(m>n) coefficient matrix(Rank Defect);

* b ---> The n-by-1 right-hand side vector;

* alpha : regularization factor;

* x ---> The n-by-1 solution vector.

template

Vector > tikhonov( const Matrix > &A,

const Vector > &b,

Type &alpha )

{

int m = A.rows(),

n = A.cols();

assert( m == b.size() );

CSVD svd;

svd.dec( A );

Matrix > U = svd.getU();

Matrix > V = svd.getV();

Vector s = svd.getSV();

Vector > x(n);

int p = s.size();

Type alpha2 = alpha*alpha;

for( int i=0; i

s[i] += alpha2/s[i];

// y = U^H * b

Vector > y(p);

for( int i=0; i

for( int j=0; j

y[i] += conj(U[j][i])*b[j];

// y = y / s

for( int i=0; i

y[i] /= s[i];

// x = V * y

for( int i=0; i

for( int j=0; j

x[i] += V[i][j]*y[j];

// for( int i=0; i

// for( int j=0; j

// {

// complex sum = 0;

// for( int k=0; k

// sum += conj(U[k][j])*b[k];

// x[i] += sum * V[i][j] / s[j];

// }

return x;

}

运行结果：

The original matrix A : size: 8 by 6

64.0000 2.0000 3.0000 61.0000 60.0000 6.0000

9.0000 55.0000 54.0000 12.0000 13.0000 51.0000

17.0000 47.0000 46.0000 20.0000 21.0000 43.0000

40.0000 26.0000 27.0000 37.0000 36.0000 30.0000

32.0000 34.0000 35.0000 29.0000 28.0000 38.0000

41.0000 23.0000 22.0000 44.0000 45.0000 19.0000

49.0000 15.0000 14.0000 52.0000 53.0000 11.0000

8.0000 58.0000 59.0000 5.0000 4.0000 62.0000

The constant vector b : size: 8 by 1

260.0000

The solution of Truncated SVD (is consistent with Matlab) :

size: 6 by 1

1.1538

1.4615

1.3846

1.4615

1.1538

The solution of Damped SVD with alpha = 0.0000 :

size: 6 by 1

1.1539

1.4615

1.3847

1.4615

1.1538

The solution of Tikhonov Regularization with alpha = 0.0000 :

size: 6 by 1

1.1538

1.4615

1.3846

1.4615

1.1538

The original complex matrix cA is: A - jA

The constant complex vector cb is: b

The solution of Truncated SVD (is consistent with Matlab) :

size: 6 by 1

(0.5769,0.5769)

(0.7308,0.7308)

(0.6923,0.6923)

(0.7308,0.7308)

(0.5769,0.5769)

The solution of Damped SVD with alpha = 0.0000 :

size: 6 by 1

(0.5769,0.5769)

(0.7308,0.7308)

(0.6923,0.6923)

(0.7308,0.7308)

(0.5769,0.5769)

The solution of Tikhonov Regularization with alpha = 0.0000 :

size: 6 by 1

(0.5769,0.5769)

(0.7308,0.7308)

(0.6923,0.6923)

(0.7308,0.7308)

(0.5769,0.5769)

Process returned 0 (0x0) execution time : 0.125 s

Press any key to continue.

weixin_39939918

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
c语言求解病态方程组,病态方程组求解方法的C++实现（支持复系数方程组）

头文件：/** Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com** This program is free software; you can redistribute it and/or modify it* under the terms of the GNU General Public License as p...
复制链接

扫一扫