高斯消元N×M的矩阵 [英] Gauss Elimination for NxM matrix
本文介绍了高斯消元N×M的矩阵的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
/ *程序来证明高斯＆LT;强烈的阶级=高亮＆GT;消除＆LT; / STRONG＆GT;
在一组线性方程组的
* /
＃包括＆LT;的iostream＆GT;
＃包括＆LT; CMATH＆GT;
＃包括＆LT;载体＆GT;
使用名字空间std;
常量双EPS = 1.E-15;
/ * preliminary旋转策略
旋转功能
* /
双枢轴（矢量＆lt;矢量＆lt;双＆GT;＆GT;＆功放;一，矢量＆lt;双＆GT;和b，int i）以
{
INT N = a.size（）;
INT J =;
双吨= 0;
为（中间体K = 1; K＆n种; K + = 1）
{
双AKI =晶圆厂（一个[k]的[I]）;
如果（AKI＆GT; T）
{
T = AKI;
J = K表;
}
}
如果（J＆GT;ⅰ）
{
双模拟;
为（中间体L = 0; l n种; L + = 1）
{
假= A [1] [1];
A [1] [1] = A [J] [L]
一个[J] [L] =假;
}
双温度= B [J]。
B〔I] = B [J]。
B〔J] =温度;
}
返回[I] [I];
}
/ *转发＆LT;强烈的阶级=高亮＆GT;消除＆LT; / STRONG＆GT; * /
无效triang（矢量＆lt;矢量＆lt;双＆GT;＆GT;＆功放;一，矢量＆lt;双＆GT;和b）
{
INT N = a.size（）;
对（INT I = 0; I＆n种-1; I + = 1）
{
双DIAG =枢轴（A，B，I）;
如果（晶圆厂（诊断）＆LT; EPS）
{
COUT＆LT;＆LT;零DET＆LT;＆LT; ENDL;
返回;
}
为（诠释J = i + 1的; J＆n种; J + = 1）
{
双MULT = A [J] [我] /诊断;
为（中间体K = 1 + 1; K＆n种; K + = 1）
{
一个[J] [K] - = MULT * A [1] [K]
}
B〔J] - = MULT * B [I]
}
}
}
/ *
积两个向量
* /
双dotProd（矢量＆lt;双＆GT;＆安培; U，矢量＆lt;双＆GT;＆安培; V，INT K1，K2 INT）
{
双总和= 0;
的for（int i = K1; I＆LT; = K2;我+ = 1）
{
总和+ = U [I] * V [I]
}
返回总和;
}
/ *
回代步骤
* /
无效backSubst（矢量＆lt;矢量＆lt;双＆GT;＆GT;＆功放;一，矢量＆lt;双＆GT;和B，矢量＆lt;双＆GT;＆安培; X）
{
INT N = a.size（）;
的for（int i = N-1; I＆GT; = 0;我 - = 1）
{
X [I] =（二[I] - dotProd（A [1]中，x，i + 1的，N-1））/ A [1] [I];
}
}
/ *
精炼高斯＆LT;强烈的阶级=高亮＆GT;消除＆LT; / STRONG＆GT;程序
* /
无效高斯（矢量＆lt;矢量＆lt;双＆GT;＆GT;＆功放;一，矢量＆lt;双＆GT;和B，矢量＆lt;双＆GT;＆安培; X）
{
triang（A，B）;
backSubst（A，B，x）的;
}
//例主体计划
诠释的main（）
{
INT N;
CIN＆GT;＆GT; N;
矢量＆lt;矢量＆lt;双＆GT; ＆GT;一个;
矢量＆lt;双＆GT; X;
矢量＆lt;双＆GT; B：
的for（int i = 0;我n种;我++）{
矢量＆lt;双＆GT;温度;
为（诠释J = 0; J＆n种; J ++）{
int无;
CIN＆GT;＆GT;没有;
temp.push_back（无）;
}
a.push_back（临时）;
b.push_back（0）;
x.push_back（0）;
}
/ *
的for（int i = 0;我n种;我++）{
int无;
CIN＆GT;＆GT;没有;
b.push_back（无）;
x.push_back（0）;
}
* /
高斯（A，B，x）的;
用于（为size_t I = 0; I＆LT; x.size（）;我++）{
COUT＆LT;＆LT; ×〔1] - ;＆其中; ENDL;
}
返回0;
}
以上高斯eleimination算法正常工作的N×N的矩阵。但我需要它的工作在N×M的矩阵。谁能帮我做到这一点？我不是很擅长数学。我得到了一些网站本code，我坚持它。
解决方案
- （可选）了解这个。请在纸上的一些例子。
- 的不要写$ C $下高斯消自己的。如果没有一定的照顾，天真高斯旋转不稳定。你有规模的线，并采取以最大的元素旋转的护理，出发点是有。请注意，这个建议仍然持有大部分线性代数算法。
- 如果您想解决方程组， LU分解，的QR分解（稳定的比吕，但速度较慢）， Cholesky分解 （在的情况下的系统是对称的）或 SVD 的（在系统中不是正方形的情况下）几乎总是更好的选择。高斯消元法是最适合不过的计算决定。
- 使用从 LAPACK 的算法，这就需要高斯消元法（如解决系统，或计算决定）的问题。真。不要滚你自己。既然你是做C ++，你可能会感兴趣犰狳这需要照顾很多东西给你。
- 如果你必须推出自己的教学起见，先看看的数字食谱，第3版的。第2版可以在网上免费的，如果你是低预算/有一个库中没有获得被发现。
- 作为一般的意见，不要code算法，你不明白。
/* Program to demonstrate gaussian <strong class="highlight">elimination</strong>
on a set of linear simultaneous equations
*/
#include <iostream>
#include <cmath>
#include <vector>
using namespace std;
const double eps = 1.e-15;
/*Preliminary pivoting strategy
Pivoting function
*/
double pivot(vector<vector<double> > &a, vector<double> &b, int i)
{
int n = a.size();
int j=i;
double t=0;
for(int k=i; k<n; k+=1)
{
double aki = fabs(a[k][i]);
if(aki>t)
{
t=aki;
j=k;
}
}
if(j>i)
{
double dummy;
for(int L=0; L<n; L+=1)
{
dummy = a[i][L];
a[i][L]= a[j][L];
a[j][L]= dummy;
}
double temp = b[j];
b[i]=b[j];
b[j]=temp;
}
return a[i][i];
}
/* Forward <strong class="highlight">elimination</strong> */
void triang(vector<vector<double> > &a, vector<double> &b)
{
int n = a.size();
for(int i=0; i<n-1; i+=1)
{
double diag = pivot(a,b,i);
if(fabs(diag)<eps)
{
cout<<"zero det"<<endl;
return;
}
for(int j=i+1; j<n; j+=1)
{
double mult = a[j][i]/diag;
for(int k = i+1; k<n; k+=1)
{
a[j][k]-=mult*a[i][k];
}
b[j]-=mult*b[i];
}
}
}
/*
DOT PRODUCT OF TWO VECTORS
*/
double dotProd(vector<double> &u, vector<double> &v, int k1,int k2)
{
double sum = 0;
for(int i = k1; i <= k2; i += 1)
{
sum += u[i] * v[i];
}
return sum;
}
/*
BACK SUBSTITUTION STEP
*/
void backSubst(vector<vector<double> > &a, vector<double> &b, vector<double> &x)
{
int n = a.size();
for(int i = n-1; i >= 0; i -= 1)
{
x[i] = (b[i] - dotProd(a[i], x, i + 1, n-1))/ a[i][i];
}
}
/*
REFINED GAUSSIAN <strong class="highlight">ELIMINATION</strong> PROCEDURE
*/
void gauss(vector<vector<double> > &a, vector<double> &b, vector<double> &x)
{
triang(a, b);
backSubst(a, b, x);
}
// EXAMPLE MAIN PROGRAM
int main()
{
int n;
cin >> n;
vector<vector<double> > a;
vector<double> x;
vector<double> b;
for (int i = 0; i < n; i++) {
vector<double> temp;
for (int j = 0; j < n; j++) {
int no;
cin >> no;
temp.push_back(no);
}
a.push_back(temp);
b.push_back(0);
x.push_back(0);
}
/*
for (int i = 0; i < n; i++) {
int no;
cin >> no;
b.push_back(no);
x.push_back(0);
}
*/
gauss(a, b, x);
for (size_t i = 0; i < x.size(); i++) {
cout << x[i] << endl;
}
return 0;
}
The above gaussian eleimination algorithm works fine on NxN matrices. But I need it to work on NxM matrix. Can anyone help me to do it? I am not very good at maths. I got this code on some website and i am stuck at it.
解决方案
- (optional) Understand this. Do some examples on paper.
- Don't write code for Gaussian elimination yourself. Without some care, the naive gauss pivoting is unstable. You have to scale the lines and take care of pivoting with the greatest element, a starting point is there. Note that this advice still holds for most linear algebra algorithms.
- If you want to solve systems of equations, LU decomposition, QR decomposition (stabler than LU, but slower), Cholesky decomposition (in the case the system is symmetric) or SVD (in the case the system is not square) are almost always better choices. Gaussian elimination is best for computing determinants however.
- Use the algorithms from LAPACK for the problems which need Gaussian elimination (eg. solving systems, or computing determinants). Really. Don't roll your own. Since you are doing C++, you may be interested in Armadillo which takes care of a lot of things for you.
- If you must roll your own for pedagogical reasons, have a look first at Numerical Recipes, version 3. Version 2 can be found online for free if you're low on budget / have no access to a library.
- As a general advice, don't code algorithms you don't understand.
这篇关于高斯消元N×M的矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
