crout.cpp

来自「解线性方程组的ＣＲＯＵＴ分解法的Ｃ语言程序实现．」· C++ 代码 · 共 44 行

#include<iostream>
using namespace std;
void main(){
double a[4][4]=
	{ {0.94455,-2.81716,-2.25806,-2.89224},
	{-6.20716,-3.85174,9.77783,3.49712},
	{-7.56157,-4.03668,9.38475,-8.64498},
	{-8.94467,-5.31114,-0.98666,-4.10932} };
double b[4]={1.51158,-6.46352,1.23264,0.61068};
 double l[4][4],u[4][4];
 double x[4],y[4];
 for(int i=0;i<4;i++){
  l[i][0]=a[i][0];
  u[0][i]=a[0][i]/l[0][0];
 }
 for(int k=1;k<4;k++){
  u[k][k]=1;
  for(int i=k;i<4;i++){
   l[i][k]=a[i][k];
   for(int r=0;r<k;r++)
    l[i][k]-=l[i][r]*u[r][k];
  }
  for(int j=k+1;j<4;j++){
   u[k][j]=a[k][j]/l[k][k];
   for(int r=0;r<k;r++)
    u[k][j]-=l[k][r]*u[r][j]/l[k][k];
  }
 }
 y[0]=b[0]/l[0][0];
 for(int m=1;m<4;m++){
  y[m]=b[m]/l[m][m];
  for(int r=0;r<m;r++)
   y[m]-=l[m][r]*y[r]/l[m][m];
 }
 x[3]=y[3];
 for(int n=2;n>=0;n--){
  x[n]=y[n];
  for(int r=n+1;r<4;r++)
   x[n]-=u[n][r]*x[r];
 }
 cout<<"方程组的解"<<endl;
 for(int g=0;g<4;g++)
  cout<<x[g]<<" "<<endl;
}