LU演算法設計
⑴ 求C語言課程設計：用高斯列主元消元法解線性方程組
這里向你推薦一下克魯特演算法（其實就是對高斯列主元消元法進行優化，使之更適合於計算機編程），首先將矩陣A進行LU分解（將系數矩陣分解成一個上三角矩陣和一個下三角矩陣），分解的過程中用到了隱式的主元尋找法，同時利用克魯特演算法可以將兩個n*n矩陣壓縮到一個n*n矩陣中，大大節省了存儲空間提高了計算速度。
方程可化為L*U*x=B，令U*x=y --->L*y=B
然後利用回代先求y，再利用y求x
因為該方法在求解過程中不涉及增廣矩陣所以矩陣B幾乎不參與什麼運算，所以它的計算速度應該能夠達到高斯列主元消元法的三倍，但原理與其基本一致。
而且我在程序中使用了動態數組方便你今後進行擴展。
以下程序按照《矩陣論第二版》和《C語言數值計演算法方法大全》編寫，LU分解部分程序主要參考了《C語言數值計演算法方法大全》第二章的程序
如果你需要詳細的理論講解我可以將這兩本書和源程序發給你，上面的論述相當詳細足夠你答辯用的了，我的郵箱[email protected]
計算結果：
A矩陣:
2 2 5
3 4 7
1 3 3
B矩陣:
5
6
5
解矩陣:
x 1=-7
x 2=0.333333
x 3=3.66667
Press any key to continue
#include <cmath>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <functional>
#include <vector>
#include <algorithm>
using namespace std;
#define TINY 1.0e-20 //A small number.
#define N 3
void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d);//對矩陣進行LU分解
void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b);//對矩陣進行前向和後向回代
void root(vector<vector<float> > &x,vector<float> &col);//解方程結果保存在y中
void iniv(vector<vector<float> > &x,vector<float> line,int n);//對二維動態數組進行初始化
void main()
{
int i,j,n=N;//輸入矩陣的維數
float A[N][N]={{2,2,5},{3,4,7},{1,3,3}};//左邊A矩陣
float B[N]={5,6,5};//右邊B矩陣
vector<vector<float> > x;//建立動態二維數組存放A，保證你的程序進行擴展時只改A，B，N
vector<float> line;
vector<float> y(n);//建立動態數組存放B
iniv(x,line,n);
y.clear();
for(i=0;i<n;i++)//將A賦給x，B賦給y
{
y.push_back(B[i]);
for(j=0;j<n;j++)
{
x[i].push_back(A[i][j]);
}
}
cout<<"A矩陣:"<<endl;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<x[i][j]<<" ";
}
cout<<endl;
}
cout<<"B矩陣:"<<endl;
for(i=0;i<n;i++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<y[i]<<endl;
}
root(x,y);//求根
cout<<"解矩陣:"<<endl;
for(i=0;i<n;i++)
{
cout<<setw(2)<<setiosflags(ios::left)<<"x"<<i+1<<"="<<y[i]<<endl;
}
cout<<endl;
}
void root(vector<vector<float> > &x,vector<float> &col)
{
int n=x.size(),i=0,j=0;
vector<int> index(n);//用於記錄尋找主元素過程中對矩陣的初等變換
index.clear();
float m=1.0;//記錄變換方式，此程序中無用
ludcmp(x,n,index,m);//進行LU分解
lubksb(x,n,index,col);//根據分解結果進行回帶
}
//以下程序按照《矩陣論第二版》和《C語言數值計演算法方法大全》編寫，LU分解部分程序主要參考了《C語言數值計演算法方法大全》第二章的程序
//如果你需要詳細的理論講解我可以將這兩本書和源程序發給你，我的郵箱[email protected]
void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d)
{
int i,imax,j,k;
float big=0,m=0,sum=0,temp=0;
vector<float> vv(n);
vv.clear();
d=1.0;
for (i=0;i<n;i++)
{
big=0.0;
for (j=0;j<n;j++)
if ((temp=fabs(a[i][j])) > big)
big=temp;
vv[i]=1.0/big;
}
for (j=0;j<n;j++)
{
for (i=0;i<j;i++)
{
sum=a[i][j];
for (k=0;k<i;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<n;i++)
{
sum=a[i][j];
for (k=0;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (m=vv[i]*fabs(sum)) >= big)
{
big=m;
imax=i;
}
}
if (j != imax)
{
for (k=0;k<n;k++)
{
m=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=m;
}
d = -(d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0)
a[j][j]=TINY;
if (j != n)
{
m=1.0/(a[j][j]);
for (i=j+1;i<n;i++)
a[i][j] *= m;
}
}
}
void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b)
{
int i,ii=0,ip,j;
float sum;
for(i=0;i<n;i++)//按LU分解時尋找主元所進行的初等變換進行反邊變換。
{
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
b[i]=sum;
}
sum=0;
for (i=1;i<n;i++)
{
sum=0;
for(j=0;j<i;j++)
{
sum+=a[i][j]*b[j];
}
b[i]=b[i]-sum;
}
b[n-1]=b[n-1]/a[n-1][n-1];
for (i=n-2;i>=0;i--)
{
sum=0;
for(j=i+1;j<n;j++)
{
sum+=a[i][j]*b[j];
}
b[i]=(b[i]-sum)/a[i][i];
}
}
void iniv(vector<vector<float> > &x,vector<float> line,int n)
{
int i,j;
for(i=0;i<n;i++)
{
x.push_back(line);
for(j=0;j<n;j++)
{
x[i].clear();
}
}
}