C语言用矩阵求解方程组
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了C语言用矩阵求解方程组相关的知识,希望对你有一定的参考价值。
yanjiandong1224@163.com
19941224
程序在这个邮箱里,行列式的计算没有问题,好像问题出在了把行列式变换上,就是slove(1,n)这个函数的运算上,那位大神可以看看问题到底出在哪…… 不胜感激!!!
//作者:与你看日出
//日期:2009年4月25日 星期六
//说明:输出值只能是小数(最多六位),如x=1.876546
//Han: 初始设的最多未知数的个数,运行程序后只能比它的个数小
//hang:计算中途中自己改变的未知数的个数
//JUZHEN:初始的示例矩阵
#include <stdio.h>
#include <math.h>
#define Han 200//(可自设)多元一次方程组有n行n+1列(多的一列是等号右边的值),给出行数就能确定矩阵,
#define JUZHEN 1,1,1,1,5,1,2,-1,4,-2,2,-3,-1,-5,-2,3,1,2,11,0//示例一个
main()
int i,j,k,m,n,t,cf,hang=4;
float temp;
float AA[Han][Han+1]=JUZHEN;//定义所要计算的数组
do//判断是否重试
for(i=0;i<hang;i++)//输出所定义的数组
printf("\\n");
for(j=0;j<hang+1;j++)
printf("%g\\t",AA[i][j]);
printf("\\n");
printf("是否自己输入?是:1;否:0");
scanf("%d",&t);
if(t==1)//判断是否自己输入数组
printf("输入未知数的个数");
scanf("%d",&hang);
for(i=0;i<hang;i++)//输入所定义的数组
for(j=0;j<hang+1;j++)
printf("第%d行第%d列的数为:",i+1,j+1);
scanf("%f",&AA[i][j]);
for(k=0;k<hang;k++)//这个大循环将数组的左下角转化为0
while(AA[k][k]==0)//如果第K行K列的那个数为0,则加和重组一行。
for(m=k+1;m<hang;m++)
for(n=k;n<hang+1;n++)
AA[k][n]+=AA[m][n];
for(i=k;i<hang;i++)//将第K列下面变为1
temp=AA[i][k];
for(j=k;j<hang+1;j++)//将每列变为1
AA[i][j]/=temp;
for(i=k+1;i<hang;i++)//将下面的数列与上面的数列相减使其下面为0
for(j=0;j<hang+1;j++)
AA[i][j]-=AA[k][j];
for(k=hang-2;k>=0;k--)//这个大循环将数组的右上角转化为0
for(i=k+1;i<hang+1-1;i++)//将第i列上面变为0
AA[k][hang+1-1]-=AA[k][i]*AA[i][hang+1-1];
AA[k][i]=0;
for(i=0;i<hang;i++)//输出该矩阵(也就是多元一次方程组)的解
printf("\\n");
for(j=0;j<hang+1;j++)
printf("%g\\t",AA[i][j]);
printf("\\n未知数的值为:\\n");
for(i=0;i<hang;i++)//输出该矩阵(也就是多元一次方程组)的解
printf("x(%d)=\\t%g\\n",i+1,AA[i][hang+1-1]);
printf("\\n");
printf("是否再试一次?是:1;否:0");
scanf("%d",&cf);
while(cf==1);//判断是否重试
追问
可以帮忙看一下我的问题在哪里吗?
参考技术A (1)由a2,6,a3成等差数列,得12=a2+a3…(2分)
又an为等比数列,且a1=2,
故12=2q+2q2…(3分)
解得q=2,或q=-3,
又q>0…(5分),
∴q=2,
∴an=2?2n?1=2n…(7分)
(2)∵bn=log22n=n,
用 GSL 求解超定方程组及矩阵的奇异值分解(SVD)
用 GSL 求解超定方程组及矩阵的奇异值分解(SVD)
最近在学习高动态图像(HDR)合成的算法,其中需要求解一个超定方程组,因此花了点时间研究了一下如何用 GSL 来解决这个问题。
GSL 里是有最小二乘法拟合(Least-Squares Fitting)的相关算法,这些算法的声明在 gsl_fit.h 中,所以直接用 GSL 提供的 gsl_fit_linear 函数就能解决这个问题。不过我想顺便多学习一些有关 SVD 的知识。所以就没直接使用 gsl_fit_linear 函数。
SVD 分解的一些基本概念
关于 SVD 有两篇不错的科普文:
- A Singularly Valuable Decomposition: The SVD of a Matrix
- We Recommend a Singular Value Decomposition
建议大家找来读读,这两篇文章似乎都已经有人翻译成中文了。
所谓 SVD,就是把一个矩阵
上面式子中的
矩阵
因为矩阵
举例来说,如果矩阵
这里
GSL 中的相关函数
gsl 中提供了好几个函数来计算 SVD:
- gsl_linalg_SV_decomp 这个是最基本的,使用 Golub-Reinsch SVD 算法,一般我们用这个就够了。
- gsl_linalg_SV_decomp_mod 这个是改进后的 Golub-Reinsch SVD 算法,当
M?N 时比 Golub-Reinsch SVD 算法要快。 - gsl_linalg_SV_decomp_jacobi 这个算法用到了 Jacobi 正交化,号称计算结果比 Golub-Reinsch SVD 算法要更准确。
除此之外,还有个 gsl_linalg_SV_solve 函数。这个就是利用 SVD 的结果来求解线性代数方程组的。
把这几个函数组合一下就可以合成一个求解线性代数方程组
下面是函数代码:
void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x)
{
int rows = A->size1;
int cols = A->size2;
gsl_vector * work = gsl_vector_alloc (cols);
gsl_vector * S = gsl_vector_alloc (cols);
gsl_matrix * U = gsl_matrix_alloc(rows, cols);;
gsl_matrix * V = gsl_matrix_alloc(cols, cols);
gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中
gsl_linalg_SV_decomp( U, V, S, work );
gsl_linalg_SV_solve ( U, V, S, b, x );
gsl_vector_free(work);
gsl_vector_free(S);
gsl_matrix_free(V);
gsl_matrix_free(U);
}
当
下面举一个具体的例子:
下面是测试代码:
void test1()
{
double a_data[] = {1.4, 2.1, 2.1, 7.4, 9.6,
1.6, 1.5, 1.1, 0.7, 5.0,
3.8, 8.0, 9.6, 5.4, 8.8,
4.6, 8.2, 8.4, 0.4, 8.0,
2.6, 2.9, 0.1, 9.9, 7.7};
gsl_matrix_view A = gsl_matrix_view_array (a_data, 5, 5);
double b_data[] = {1.1, 1.6, 4.7, 9.1, 0.1};
gsl_vector_view b = gsl_vector_view_array (b_data, 5);
gsl_vector * x = gsl_vector_alloc (5);
linearSolve_SVD(&A.matrix, &b.vector, x);
gsl_vector_fprintf (stdout, x, "%f");
qDebug() << "";
gsl_vector * bb = gsl_vector_alloc (5);
gsl_blas_dgemv (CblasNoTrans, 1, &A.matrix, x, 0, bb);
gsl_vector_fprintf (stdout, bb, "%f");
}
输出结果如下:
-5.208566
5.736694
-2.537472
-1.029814
0.968151
1.100000
1.600000
4.700000
9.100000
0.100000
可以看出计算结果还是很准确的。
当
测试代码如下:
void test3()
{
double a_data[] = {2, 4,
3, -5,
1, 2};
gsl_matrix_view A = gsl_matrix_view_array (a_data, 3, 2);
double b_data[] = {11, 3, 6};
gsl_vector_view b = gsl_vector_view_array (b_data, 3);
gsl_vector * x = gsl_vector_alloc (2);
linearSolve_SVD(&A.matrix, &b.vector, x);
gsl_vector_fprintf (stdout, x, "%f");
qDebug() << "";
gsl_vector * bb = gsl_vector_alloc (3);
gsl_blas_dgemv (CblasNoTrans, 1, &A.matrix, x, 0, bb);
gsl_vector_fprintf (stdout, bb, "%f");
}
计算结果如下:
3.090909
1.254545
11.200000
3.000000
5.600000
如果
方程很简单,口算就可以出结果,这个方程的解是:
下面用我们的代码计算一下。
void test4()
{
double a_data[] = {1, 2,
2, 4};
gsl_matrix_view A = gsl_matrix_view_array (a_data, 2, 2);
double b_data[] = {3, 6};
gsl_vector_view b = gsl_vector_view_array (b_data, 2);
gsl_vector * x = gsl_vector_alloc (2);
linearSolve_SVD(&A.matrix, &b.vector, x);
gsl_vector_fprintf (stdout, x, "%f");
qDebug() << "";
gsl_vector * bb = gsl_vector_alloc (2);
gsl_blas_dgemv (CblasNoTrans, 1, &A.matrix, x, 0, bb);
gsl_vector_fprintf (stdout, bb, "%f");
}
结果是:
-3.400000
3.200000
3.000000
6.000000
可以验算,
下面我们将 SVD 相关的功能封装成一个类,以方便我们提取
另外,当我们一个
头文件如下:
#ifndef GSLSINGULARVALUEDECOMPOSITION_H
#define GSLSINGULARVALUEDECOMPOSITION_H
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>
void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x);
class GslSVD
{
public:
GslSVD();
~GslSVD();
int SV_decomp(const gsl_matrix * A);
int SV_decomp_mod(const gsl_matrix * A);
int SV_decomp_jacobi (gsl_matrix * A);
int SV_solve(const gsl_vector *b, gsl_vector *x);
gsl_vector * getVectorS();
gsl_matrix * getMatrixU();
gsl_matrix * getMatrixV();
int trimVectorS(double abseps);
private:
gsl_vector * S;
gsl_matrix * U;
gsl_matrix * V;
void alloc_suv(int rows, int cols);
};
#endif // GSLSINGULARVALUEDECOMPOSITION_H
cpp 文件如下:
#include "gsl_SVD.h"
void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x)
{
int rows = A->size1;
int cols = A->size2;
gsl_vector * work = gsl_vector_alloc (cols);
gsl_vector * S = gsl_vector_alloc (cols);
gsl_matrix * U = gsl_matrix_alloc(rows, cols);;
gsl_matrix * V = gsl_matrix_alloc(cols, cols);
gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中
gsl_linalg_SV_decomp( U, V, S, work );
gsl_linalg_SV_solve ( U, V, S, b, x );
gsl_vector_free(work);
gsl_vector_free(S);
gsl_matrix_free(V);
gsl_matrix_free(U);
}
int GslSVD::trimVectorS(double abseps)
{
int count = 0;
for(int i = 0; i < S->size; i++)
{
if(fabs(gsl_vector_get(S, i)) < abseps)
{
count ++;
gsl_vector_set(S, i, 0);
}
}
return count;
}
gsl_vector * GslSVD::getVectorS()
{
if(S == NULL) return NULL;
gsl_vector * s = gsl_vector_alloc(S->size);
gsl_vector_memcpy(s, S);
return s;
}
gsl_matrix * GslSVD::getMatrixU()
{
if(U == NULL) return NULL;
gsl_matrix * u = gsl_matrix_alloc(U->size1, U->size2);
gsl_matrix_memcpy(u, U);
return u;
}
gsl_matrix * GslSVD::getMatrixV()
{
if(V == NULL) return NULL;
gsl_matrix * v = gsl_matrix_alloc(V->size1, V->size2);
gsl_matrix_memcpy(v, V);
return v;
}
GslSVD::GslSVD()
{
S = NULL;
U = NULL;
V = NULL;
}
void GslSVD::alloc_suv(int rows, int cols)
{
if( S != NULL )
{
gsl_vector_free(S);
gsl_matrix_free(U);
gsl_matrix_free(V);
}
S = gsl_vector_alloc (cols);
U = gsl_matrix_alloc(rows, cols);
V = gsl_matrix_alloc(cols, cols);
}
int GslSVD::SV_decomp(const gsl_matrix * A)
{
int rows = A->size1;
int cols = A->size2;
gsl_vector * work = gsl_vector_alloc (cols);
alloc_suv(rows, cols);
gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中
int ret = gsl_linalg_SV_decomp( U, V, S, work );
gsl_vector_free(work);
return ret;
}
int GslSVD::SV_decomp_mod(const gsl_matrix * A)
{
int rows = A->size1;
int cols = A->size2;
gsl_vector * work = gsl_vector_alloc (cols);
gsl_matrix *X = gsl_matrix_alloc(cols, cols);
alloc_suv(rows, cols);
gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中
int ret = gsl_linalg_SV_decomp_mod( U, X, V, S, work );
gsl_matrix_free(X);
gsl_vector_free(work);
return ret;
}
int GslSVD::SV_decomp_jacobi (gsl_matrix * A)
{
int rows = A->size1;
int cols = A->size2;
alloc_suv(rows, cols);
gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中
int ret = gsl_linalg_SV_decomp_jacobi( U, V, S );
return ret;
}
int GslSVD::SV_solve(const gsl_vector *b, gsl_vector *x)
{
if(U != NULL)
{
return gsl_linalg_SV_solve (U, V, S, b, x);
}
return -1;
}
GslSVD::~GslSVD()
{
if(S != NULL)
{
gsl_vector_free(S);
gsl_matrix_free(V);
gsl_matrix_free(U);
}
}
下面用这个类来计算一下刚才的问题:
void test5()
{
double a_data[] = {1, 2,
2, 4};
gsl_matrix_view A = gsl_matrix_view_array (a_data, 2, 2);
GslSVD svd;
svd.SV_decomp(&A.matrix);
puts("S = ");
gsl_vector_fprintf (stdout, svd.getVectorS(), "%f");
puts("\nV = ");
gsl_matrix_fprintf (stdout, svd.getMatrixV(), "%f");
double b_data[] = {3, 6};
gsl_vector_view b = gsl_vector_view_array (b_data, 2);
gsl_vector * x = gsl_vector_alloc (2);
svd.SV_solve(b, x);
puts("\nx = ");
gsl_vector_fprintf (stdout, x, "%f");
}
结果如下:
S =
5.000000
0.000000
V =
-0.447214
-0.894427
-0.894427
0.447214
x =
-3.400000
3.200000
我们注意到
大家可以验证一下,这个解是正确的。
另外,我写的类中还提供了一个 trimVectorS(double abseps) 函数,利用这个函数,可以将
好了,就先写这么多。希望对大家有用。
以上是关于C语言用矩阵求解方程组的主要内容,如果未能解决你的问题,请参考以下文章