线性方程组AX=B的数值计算方法实验

学号：姓名：梁哲豪

一、 实验描述

在自然科学和工程技术中很多问题的解决常常归结为解线性代数方程组。例如电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组问题，用差分法或者有限元法解常微分方程，偏微分方程边值问题等都导致求解线性方程组，而且后面几种情况常常归结为求解大型线性方程组。线性代数方面的计算方法就是研究求解线性方程组的一些数值解法与研究计算矩阵的特征值及特征向量的数值方法。 关于线性方程组的数值解法一般有两类：

直接法：若在计算过程中没有舍入误差，经过有限步算术运算，可求得方程组的精确解的方法。

迭代法：用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有占存储单元少，程序设计简单，原始系数矩阵在迭代过程中不变等优点，但存在收敛性及收敛速度等问题。 上三角线性方程组的求解： 基本算法：

高斯消元法：将原方程组化为三角形方阵的方程组：

（k=1,2,…,n-1; i=k+1,k+2, …,n ;j=k+1,k+2, …,n+1） 由回代过程求得原方程组的解：

LU分解法：

将系数矩阵A转化为A=L*U，L为单位下三角矩阵，U为普通上三角矩阵，然后通过解方程组l*y=b，u*x=y,来求解x。

二、实验内容

1、许多科学应用包含的矩阵带有很多零。在实际情况中很重要的三角形线性方程组有如下形式： ……

构造一个程序求解三角形线性方程组。可假定不需要变换。而且可用第k行消去第k+1行的xk。 核心代码：

#include #include #include #define N 4//矩阵阶数

void ColPivot(double c[N][N+1],double []); //函数声明 void main(){ int i,j; double x[N];

double c[N][N+1]={1,3,5,7,1, 2,-1,3,5,2, 0,0,2,5,3, -2,-6,-3,1,4};

cout<<\"----------------------------------------\"<cout<cout<<\"右侧矩阵 y 为: \\n\"; for(i=0;icout<cout<<\"----------------------------------------\"<cout<<\"变换后得到的三角矩阵: \\n\"; for(i=0;icout<cout<<\"变换后的右侧矩阵 y 为: \\n\"; for(i=0;icout<cout<<\"----------------------------------------\"<cout<<\" x[\"<cout<<\"----------------------------------------\"<void ColPivot(double c[N][N+1],double x[]) {

int i,j,k; double p,max; double t[N];

for(i=0;i<=N-2;i++) { max=0; k=i;

for(j=i+1;jmax) { k=j;

max=fabs(c[j][i]); //选主元 }

if(k!=i)

for(j=i;j<=N;j++) {

p=c[i][j];

c[i][j]=c[k][j]; //选出主元后进行交换 c[k][j]=p; }

for(j=i+1;jp=c[j][i]/c[i][i]; for(k=i;k<=N;k++)

c[j][k]-=p*c[i][k]; //高斯消去，进行计算 }

}

for(i=0;ifor(i=N-1;i>=0;i--) //利用回代法求最终解 {

for(j=N-1;j>i;j--) t[i]-=c[i][j]*x[j]; x[i]=t[i]/c[i][i]; } }

运行结果：（以具体方程组为例）

2、（PA=LU：带选主元的分解法）求解线性方程组AX=B,其中：

A= B=

核心代码：

#include #include #define L 30

double a[L][L],b[L],l[L][L],u[L][L],x[L],y[L]; int main() {

int n,i,j,k,r;

printf(\"请输入矩阵元次：\\n\"); scanf(\"%d\",&n);

printf(\"请输入矩阵各项：\\n\"); for(i=1;i<=n;++i) {

for(j=1;j<=n;++j) {

scanf(\"%lf\",&a[i][j]); } }

printf(\"请输入方程组的常数项：\\n\"); for(i=1;i<=n;++i) {

scanf(\"%lf\",&b[i]); }

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

for(j=1;j<=n;++j) {

l[i][j]=0; u[i][j]=0.0; } }

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

for(j=k;j<=n;++j) {

u[k][j]=a[k][j]; for(r=1;ru[k][j]-=l[k][r]*u[r][j]; } }

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

l[i][k]=a[i][k]; for(r=1;rl[i][k]-=l[i][r]*u[r][k]; }

l[i][k]/= u[k][k]; }

l[k][k]=1.0; }

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

y[i] = b[i]; for(j=1;jy[i]-=l[i][j]*y[j]; } }

for(i=n;i>0;--i) { x[i] = y[i];

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

x[i]-=u[i][j]*x[j]; }

x[i]/= u[i][i]; }

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

printf(\"%0.2lf\\n\",x[i]); } return 0; }

运行结果：

3、使用程序3.3求解线性方程组AX=B，其中，A= [aij] N×N= i j-1,而且B=[b ij] N×1, b11=N,当i≥2时，bi1=(iN-1)/(i-1),对N=3,7,11的情况分别求解。精确解为X=[1 1 1 …1 1]’。对得到的结果与精确解的差异进行解释。 核心代码：

function X=lufact(A,B)

%input--A is an N*N nonsingular matrix % -E is an N*N eyes

%output-m is an N*N invmatrix of A

%initialize X, Y, the temporary storage matrix C, and the row %permutation information matrix t [N,N]=size(A); X=zeros(N,1); Y=zeros(N,1); C=zeros(1,N); R=1:N; for p=1:N-1

%find the privot row for column p [max1,j]=max(abs(A(p:N,p))); %interchange row p and j C=A(p,:);

A(p,:)=A(j+p-1,:); A(j+p-1,:)=C; d=R(p); R(p)=R(j+p-1); R(j+p-1)=d; if A(p,p)==0

'A is singular. No unique solution' break end

%calculate multiplier and place in subdiagonal protion of A for k=p+1:N

mult=A(k,p)/A(p,p); A(k,p)=mult;

A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N); end end

%solve for Y Y(1)=B(R(1)); for k=2:N

Y(k)=B(R(k))-A(k,1:k-1)*Y(1:k-1); end

%solve for X X(N)=Y(N)/A(N,N); for k=N-1:-1:1

X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k); End

4、修改程序3.3，使得它可以通过重复求解N个线性方程组

ACJEJ 其中J=1，2，…，N

来得到A1，

则 AC1C2...CNE1E2...EN、

而且 A1C1C2...CN，保证对LU分解只计算一次。 修改程序3.3 主要变动为列矩阵B变成了单位矩阵E 算法：

（1） 输入A,E;[N,N]=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);t=1:N;

p=1。

（2） 判断p=1>N-1是否成立，若成立，执行步骤（6）；若不成立，执行步骤（3）。 （3） [max1,j]=max(abs(A(p:N,p)));C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=t(p

);t(p)=t(j+p-1);t(j+p-1)=d; （4） 判断A(p,p)==0是否成立，若成立，执行步骤（12）；若不成立，k=p+1，

执行步骤（5）。

（5） 判断k>N是否成立，若成立，返回到步骤（2）；若不成立，mult=A(k,p)/A(p,p);

A(k,p)=mult;A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);k=k+1;返回到步骤（5）。 （6） J=1。

（7） 判断J>N是否成立，若成立，执行步骤（12）；若不成立，k=2, Y(1)=E(t(1),J);

执行步骤（8）

（8） 判断k>N是否成立，若成立，执行步骤（9）；若不成立，

Y(k)=E(t(k),J)-A(k,1:k-1)*Y(1:k-1)，k=k+1,返回步骤（8）。 （9） X(N)=Y(N)/A(N,N); k=N-1.

（10）判断k> 1是否成立，若成立，执行步骤（11）；若不成立，

X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);

k=k-1;执行步骤（10）.

（11）m(1:N,J)=X(1:N);返回步骤（7）. （12）输出结果。 核心代码：

function m=newlufact(A,E)

%input--A is an N*N nonsingular matrix % -E is an N*N eyes

%output-m is an N*N invmatrix of A

%initialize X, Y, the temporary storage matrix C, and the row %permutation information matrix t [N,N]=size(A); X=zeros(N,1); Y=zeros(N,1); C=zeros(1,N); t=1:N; for p=1:N-1

%find the privot row for column p [max1,j]=max(abs(A(p:N,p))); %interchange row p and j C=A(p,:);

A(p,:)=A(j+p-1,:); A(j+p-1,:)=C; d=t(p); t(p)=t(j+p-1); t(j+p-1)=d; if A(p,p)==0

'A is singular. No unique solution' break end

%calculate multiplier and place in subdiagonal protion of A for k=p+1:N

mult=A(k,p)/A(p,p); A(k,p)=mult;

A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N); end end

%solve for Y

for J=1:N Y(1)=E(t(1),J); for k=2:N

Y(k)=E(t(k),J)-A(k,1:k-1)*Y(1:k-1); end

%solve for X X(N)=Y(N)/A(N,N); for k=N-1:-1:1

X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k); end

m(1:N,J)=X(1:N); end