数值与符号计算__LU和QR分解实验报告

数值与符号计算 实验（丑）

线性代数方程组求解

姓名： 学号： 学院：

1 实验题目

用C/C++语言实现如下函数：

bool lu(double* a, int* pivot, int n)；矩阵的LU分解

bool guass(double const* lu, int const* p, double* b, int n);求线性代数方程组的解

void qr(double* a, double* d, int n);矩阵的QR分解

bool hshld(double const*qr, double const*d, double*b, int n); 求线性代数方程组的解

2 实验要求

撰写详细的实验报告 不必修改函数界面

用高斯列选主元消去法和矩阵QR分解两种方法求解下面3个方程组，比较这两种方法的误差。

2.1 求解方程

??1?1??2?1?3?1??4?1?5?1??6

121314151617

131415161718

141516171819

11?56?11?

?67?11?78?11?

?89?11?910?11?

?

1011?

11??1???????26????x1??1?1?????1??x??237?21????11

???????x3??348? =???111?x4??????????

9??x5??45

1????11

????????x6???5610??111?????????

11??67

??????????????????

1??01?01?

? ?11???11??0

令n?n矩阵

?1000

?

??1100??1?110Gn??

?

??1?1?1?1?

??1?1?1?1?

令n维向量

Tbn??2,1,0,?1,?2,???,4?n,2?n?

求解方程

G20x?b20

注：上面两个方程准确解都是

Tx??1,1,???,1,1?

令

???2,3,4,???,n,n?T bn

求解方程

? G20x?b20

3 实验算法结果

3.1 矩阵的Lu分解算法代码如下：

bool lu(double* a, int* pivot, int n)//矩阵LU分解 {

int i,j,k;

double max,temp;

max = 0;

temp = 0;

for (i=0; i<n-1; i++)//依次对第i列进行处理 {

// 选出i列的主元,记录主元位置

max = fabs(a[n*i + i]);

pivot[i]=i;

for(j=i+1; j<n; j++)//j表示第j行

{

if( fabs(a[n*j + i])>max)

{

max= fabs(a[n*j + i]) ;

pivot[i]=j;

}

}

// 对第i列进行行变换，使得主元在对角线上 if(pivot[i]!=i)

{

for(j=i; j<n; j++)//ij与pivot[i]j换 只用对上三角进行处理 {

temp=a[n*i + j];

a[n*i + j]=a[n*pivot[i]+ j];

a[n*pivot[i]+ j]=temp;

}

}

for(j=i+1; j<n; j++)//Pi 部分下三角L

a[n*j + i]=a[n*j+i]/a[n*i+i];

for(j=i+1; j<n; j++)//计算上三角U

for(k=i+1; k<n; k++)

a[n*j + k]=a[n*j+k]- a[n*j+i]*a[n*i+k];

}

//计算下三角 L

for(i=0; i<n-2; i++)//i行k列

for(k=i+1; k<n-1;k++)

{

temp=a[n*pivot[k] + i];

a[n*pivot[k] + i]=a[k*n + i];

a[k*n + i]=temp;

}

return false ;

}

3.2 高斯求解线性方程组代码如下：

bool guass(double const* lu, int const* p, double* b, int n)//求线性代数方程组的解

{

int i,j;

double temp;

//按qivot对b行变换，与LU匹配

for(i=0; i<n-1; i++) //貌似错误在这里哦

{

temp = b[p[i]];

b[p[i]] = b[i];

b[i]=temp;

}

//Ly=b,将y的内容放入b

for(i=0; i<n; i++)

for(j=0; j<i; j++)

b[i]=b[i]-lu[n*i+j]*b[j];

//Uy=x,将x的内容放入b

for(i=n-1; i>=0; i--)

{

for(j=n-1; j>i; j--)

b[i]=b[i]-lu[n*i+j]*b[j];

b[i]=b[i]/lu[n*i+i];

}

return false;

}

3.3 矩阵的QR分解算法代码如下：

void qr(double* a, double* d, int n) //矩阵的QR分解

{

int i,j,l,k;

double tem,m;

double *temp;

temp = (double *)malloc(sizeof(double)*n);

for (i=0; i<n-1; i++)//依次对第i列进行处理

{

//获得tem值

m = 0 ;

for(j=i; j<n; j++)//j表示第j行

m = m +a[n*j + i]*a[n*j + i ];

if(a[n*i + i ]>0)

m = -sqrt(m);

else

m = sqrt(m);

//获得temp放入矩阵，并存主元d

tem = 0 ;

d[i] = m ;

a[n*i +i] = a[n*i +i] - m;

for(j=i; j<=n-1; j++)

tem=tem + a[n*j +i]*a[n*j +i];

tem= sqrt(tem);

for(j=i; j<=n-1; j++)

a[n*j +i] = a[n*j +i]/tem ;

// 调整矩阵

for(k=i+1;k<n;k++)

{

for(j=i; j<n; j++)

{

tem = 0 ;

for(l=i; l<n; l++)

tem =tem + a[n*j + i]*a[n*l + i]*a[n*l + k]; temp[j] = a[j*n+k] - 2*tem;

}

for(j=i; j<n; j++)

a[j*n+k] = temp[j];

}

}

d[n-1] = a[(n-1)*n+n-1];

}

3.4 QR求解线性方程组的代码如下：

bool householder(double const*qr, double const*d, double*b, int n)//求线性代数方程组的解

{

int i,j,l;

double rem;

double *temp;

temp = (double *)malloc(sizeof(double)*n);

for(i=0; i<n-1; i++)

{

for(j=i; j<n; j++)

{

rem = 0;

for(l=i;l<n; l++)

rem = rem + qr[l*n+i]*qr[j*n+i]*b[l];

temp[j] = b[j] - 2*rem;

}

for(j=i; j<n; j++)

b[j] = temp[j];

}

for(j=n-1; j>-1; j--)

{

for(l=n-1; l!=j;--l)

b[j] =b[j] - b[l]*qr[j*n+l];

b[j] = b[j] /d[j];

}

return false;

}