《数值分析》课程设计

问题的提出

3.3 用SOR 方法解下列方程组（去松弛因子w=1.2），要求14||||10k k X X +-∞-<。 1214

2145x x x x +=??

-=? 3.4 设 411011A ??= ???

，

计算()cond A ∞。

3.5 用选列主元Gauss 消元法求解方程组

1231231

2334721320x x x x x x x x x -+=??-+-=-??-+=?

3.6 用追赶法解三对角方程组

12345210001121000012100001210000120x x x x x -?????? ? ? ?-- ? ? ?

? ? ?=-- ? ? ?-- ? ? ? ? ? ?-????

?? 3.7 用三角分解法解方程组

123248541816862207x x x -?????? ??? ?--= ??? ? ??? ?--??????

3.8 用选主元消元法计算下列行列式

126

324951

。

一、问题分析

1. 超松弛法是迭代方法的一种加速方法，其计算公式简单，但需要选择合适的松弛因子，以保

证迭代过程有较快的收敛速度。

2. A 的条件数计算首先要获得A 的逆，而求A 的逆可以转化为求n 个方程组。

3. 完全主元消元法在计算过程中花费了大量的时间用于寻找主元。同时，各变量的位置在消元

过程中也可能会发生变化。而列选主元法则可消除这个弊病。

4. 追赶法主要是解三对角方程组。所谓追指消元过程，赶指回代过程。

5. Gauss 消元法是通过逐步消元过程，将方程组的系数矩阵A 转变为一个上三角矩阵。三角分

解法，就是把系数矩阵分解为两个三角阵。

6.将某一向量坐标同乘以某非零实数，加到另一向量上，行列式的值不变。用选主元法将行列

式矩阵变为三角阵，对角线上的数值相乘即为行列式的值。

二、编程解决

3.3Sor法c语言编程：

#include

#include

#include

#define omega 1.2 //取值不合适结果可能发散

void main()

{

double a[5][5];

double b[5],x[5],f,t,y[5]={0,0,0,0,0};

int i,j,n,cnt=0;

printf("阶数:");

scanf("%d",&n);

printf("请输入%d阶的A矩阵\n",n);

for(i=0;i

for(j=0;j

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

printf("请输入B矩阵\n");

for(i=0;i

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

printf("count\t");

for(i=0;i

printf("x[%d]\t\t",i);

printf("收敛程度\n");

do

{

for(i=0;i

x[i]=y[i];

for(i=0;i

{

t=0;

for(j=0;j

t=t+a[i][j]*(j

y[i]=x[i]+omega*(b[i]-t)/a[i][i];

printf("%d",cnt++);

for(i=0;i

printf("\t%lf",x[i]);

f=0;

for(i=0;i

f+=fabs(y[i]-x[i]);

printf("\t%g\n",f);

}while(f>1e-4 && cnt<100);

}

所得结果：

3.4 求逆、算条件数编程：

#include

#include

#include

#define N 5 //可修改，以改变可解决的最大维数。但会浪费更多的空间double a[N][N],b[N][N],c[N][N];

void input()

{

int i,j;

do

{

printf("阶数:");

scanf("%d",&n);

if(n>N || n<=0) //宏定义中N规定了数组的最大维数，超过则出错printf("溢出，请重新输入!\n");

}while(n>N || n<=0);

printf("请输入%d阶的A矩阵\n",n);

for(i=0;i

for(j=0;j

{

if(i==j)

b[i][j]=1;

else

b[i][j]=0;

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

c[i][j]=a[i][j];

}

}

void watch()

{

int i,j;

for(i=0;i

{

for(j=0;j

printf("%4.2lf\t",a[i][j]);

for(j=0;j

printf("%4.2lf\t",b[i][j]);

printf("\n");

}

system("pause");

}

void swap(double *p,double *q)

{

double tmp;

tmp=*p;

*p=*q;

}

void swap_x(int s,int t) //S<==>T

{

int j;

for(j=0;j

{

swap(&a[s][j],&a[t][j]);

swap(&b[s][j],&b[t][j]);

}

}

void self_x(int i,double s) //I=I/s

{

int j;

for(j=0;j

{

a[i][j]=a[i][j]/s;

b[i][j]=b[i][j]/s;

}

}

void inter_x(int i,int r,double t) //I=I-R*t

{

int j;

for(j=0;j

{

a[i][j]=a[i][j]-a[r][j]*t;

b[i][j]=b[i][j]-b[r][j]*t;

}

}

void sort(int k)

{

int i,maxi=k;

double max;

max=fabs(a[k][k]); //进入循环前先假定首数最大for(i=k;imax)

{

maxi=i;

max=fabs(a[i][k]);

}

/*将最大的一个数放到首数位置*/

if(maxi!=k)

swap_x(maxi,k);

}

void inv()

{

int i,k;

for(k=0;k

{

sort(k); //消元之前先排序

if(a[k][k]==0)

{

flag=0;

return;

}

self_x(k,a[k][k]);

for(i=k+1;i

inter_x(i,k,a[i][k]);

watch();

}

for(k=n-1;k>0;k--) //向上消元过程

{

for(i=0;i

inter_x(i,k,a[i][k]);

watch();

}

}

void output()

{

int i,j;

if(flag)

{

printf("A的逆矩阵如下:\n");

for(i=0;i

{

for(j=0;j

printf("%4.2lf\t",b[i][j]);

printf("\n");

}

}

else

printf("A矩阵不可逆\n");

}

double cond() //求条件数

{

double u[N],v[N],t,max1,max2;

int i,j;

for(i=0;i

{

t=0;

for(j=0;j

{

t+=fabs(b[i][j]);

}

u[i]=t;

}

for(i=0;i

{

t=0;

for(j=0;j

{

t+=fabs(c[i][j]);

}

v[i]=t;

}

max1=u[0];

max2=v[0];

for(i=0;i

{

if(u[i]>max1)

max1=u[i];

if(v[i]>max2)

max2=v[i];

}

return(max1*max2);

}

void main()

{ double cond();

input();

inv(); //求逆矩阵

output();

printf("the cond is:%lf",cond()); //求条件数printf("\n");

}

输出结果：

3.5 选列主元gauss消元法编程：#include

#include

#include

#define N 5

double a[N][N];

double b[N],x[N];

int n,flag=1;

void input()

{

int i,j;

do

{

printf("阶数:");

scanf("%d",&n);

if(n>N)

printf("溢出，请重新输入!\n");

}while(n>N);

printf("请输入%d阶的A矩阵\n",n);

for(i=0;i

for(j=0;j

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

printf("请输入B矩阵\n");

for(i=0;i

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

}

void swap(double *p,double *q)

{

double tmp;

tmp=*p;

*p=*q;

*q=tmp;

}

void swap_x(int s,int t)

{

int i;

for (i=t;i

swap(&a[s][i],&a[t][i]);

swap(&b[s],&b[t]);

}

void sort(int k)

{

int i,maxi=k;

double max;

max=fabs(a[k][k]); //进入循环前先假定首数最大for(i=k;i

if(fabs(a[i][k])>max)

{

maxi=i;

max=fabs(a[i][k]);

}

/*将最大的一个数放到首数位置*/

if(maxi!=k)

swap_x(maxi,k);

}

void gauss()

{

int i,j,k;

double t;

for(k=0;k

{

sort(k); //消元之前先排序

if(a[k][k]==0) //若发现排序后首数为0，终止!

{

flag=0;

break;

}

for(i=k+1;i

{

t=a[i][k]/a[k][k];

a[i][k]=0;

for(j=k+1;j

a[i][j]=a[i][j]-a[k][j]*t;

b[i]=b[i]-b[k]*t;

}

for(i=0;i

{

for(j=0;j

printf("%6.3lf\t",a[i][j]);

printf("%6.3lf\n",b[i]);

}

system("pause");

}

if(a[k][k]==0) //回代之前发现系数为0，终止!

flag=0;

else

{

x[n-1]=b[n-1]/a[n-1][n-1]; //回代

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

{

t=0;

for(j=i+1;j

t+=a[i][j]*x[j];

x[i]=(b[i]-t)/a[i][i];

}

}

}

void output()

{

int i;

if(flag)

for(i=0;i

printf("x[%d]=%6.3lf\t",i+1,x[i]);

else

printf("消元过程被强行终止!");

printf("\n");

}

void main()

{

input();

gauss();

output();

}

输出结果：

3.6追赶法编程：

#include

#include

#define N 5

int n;

double a[N],b[N],c[N],d[N],u[N],q[N],x[N];

void input()

{

int i;

do

{

printf("阶数:");

scanf("%d",&n);

if(n>N)

printf("溢出，请重新输入!\n");

}while(n>N);

printf("a=");

for(i=1;i

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

printf("b=");

for(i=0;i

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

printf("c=");

for(i=0;i

scanf("%lf",&c[i]);

printf("d=");

for(i=0;i

scanf("%lf",&d[i]);

}

void chase()

{

int i,j;

//变换成上三角矩阵

u[0]=c[0]/b[0];

q[0]=d[0]/b[0];

for(i=1;i

{

u[i]=c[i]/(b[i]-u[i-1]*a[i]);

q[i]=(d[i]-q[i-1]*a[i])/(b[i]-u[i-1]*a[i]);

}

for(i=0;i

for(j=0;j

{

if(ji+1)

printf("0\t");

else if(i==j)

printf("1\t");

else

printf("%4.2lf\t",u[i]);

}

printf("%lf\n",q[i]);

}

//回代过程

x[n-1]=q[n-1];

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

x[i]=q[i]-u[i]*x[i+1];

}

void output()

{

int i;

for(i=0;i

printf("x[%d]=%6.3lf\t",i,x[i]); }

void main()

{

input();

chase();

output();

}

输出结果：

3.7三角分解法解方程组

程序：

#include

#include

#define N 6 //定义一个维数N，以应用于6维以下double a[N][N],b[N],X[N],Y[N],u[N][N],v[N][N],c[N][N];

int n;

void input() //输入a系数矩阵，b常数矩阵

{

int i,j;

printf("please input the dimension(less than 7):\n"); //输入维数scanf("%d",&n);

printf("please input the matrix a[n][n]:\n");

for(i=0;i

for(j=0;j

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

printf("please input the matrix b[n]:\n");

for(i=0;i

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

}

void solveUandV() //求解u和v三角矩阵

{

int i,j,k;

double temp1,temp2;

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

{

for(j=i;j>=0;j--)

{

temp1=0;

for(k=i+1;k

{

temp1+=u[j][k]*v[k][i];

}

u[j][i]=a[j][i]-temp1;

temp2=0;

for(k=i+1;k

{

temp2+=u[i][k]*v[k][j];

}

v[i][j]=(a[i][j]-temp2)/u[i][i];

}

}

printf("上三角u:\n");

for(i=0;i

{

for(j=0;j

printf("%lf\t",u[i][j]);

printf("\n");

}

printf("下三角v:\n");

for(i=0;i

{

for(j=0;j

printf("%lf\t",v[i][j]);

printf("\n");

}

}

void substituting() //回代

{

int i,j;

double temp1,temp2;

printf("the 中间解Y is:\n");

for(i=n-1;i>=0;i--) //回代求出Y矩阵的值

{

temp1=0;

for(j=i+1;j

temp1+=Y[j]*u[i][j];

Y[i]=(b[i]-temp1)/u[i][i];

printf("%lf\t\t",Y[i]);

}

printf("\n");

for(i=0;i

temp2=0;

for(j=0;j

{

temp2+=X[j]*v[i][j];

}

X[i]=Y[i]-temp2;

}

}

void verification() //验证拆分的正确性

{

int i,j,k;

double t;

for(i=0;i

{

for(j=0;j

{

t=0;

for(k=0;k

t+=u[i][k]*v[k][j];

c[i][j]=t;

}

}

printf("the 验证的a1 is:\n"); //输出验证得到的矩阵for(i=0;i

{

for(j=0;j

printf("%lf\t",c[i][j]);

printf("\n");

}

}

void output() //输出解x

{

int i;

printf("the value of these equations is:\n");

for(i=0;i

printf("x[%d]= %lf\n",i,X[i]);

}

void main()

{

input(); //输入

solveUandV(); //拆分a为u、v上三角和下三角矩阵

substituting(); //回代求出y（中间解）和x

verification(); //验证拆分的正确性

output(); //输出解

}

结果：

3.8选主元消元法编程：

#include

#include

#include

#define N 5 //可修改，以改变可解决的最大维数。但会浪费更多的空间double a[N][N];

int n,flag=1;

void input()

{

int i,j;

do

{

printf("阶数:");

scanf("%d",&n);

if(n>N) //宏定义中N规定了数组的最大维数，超过则出错printf("溢出，请重新输入!\n");

}while(n>N);

printf("请输入%d阶的A矩阵\n",n);

for(i=0;i

for(j=0;j

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

}

void swap(double *p,double *q)

{

double tmp;

tmp=*p;

*p=*q;

*q=tmp;

}

void swap_x(int s,int t)

{

int i;

for (i=t;i

swap(&a[s][i],&a[t][i]);

}

void swap_y(int s,int t)

{

int i,tmp;

for(i=0;i

swap(&a[i][s],&a[i][t]);

}

void sort(int k)

{

int i,j,maxi=k,maxj=k;

double max;

max=fabs(a[k][k]); //进入循环前先假定首数最大for(i=k;i

for(j=k;j

if(fabs(a[i][j])>max)

{

maxi=i;

maxj=j;

max=fabs(a[i][j]);

}

/*将最大的一个数放到首数位置*/

if(maxi!=k)

swap_x(maxi,k);

if(maxj!=k)

swap_y(maxj,k);

}

void gauss()

{

int i,j,k;

double t;

for(k=0;k

{

sort(k); //消元之前先排序

if(a[k][k]==0) //若发现排序后首数为0，终止!

{

flag=0;

break;

}

for(i=k+1;i

{

t=a[i][k]/a[k][k];

a[i][k]=0;

for(j=k+1;j

a[i][j]=a[i][j]-a[k][j]*t;

}

for(i=0;i

{

for(j=0;j

printf("%6.3lf\t",a[i][j]);

printf("\n");

}

system("pause");

}

}

void determinant() //求行列式

{

double t=1.0;

int i;

for(i=0;i

{

t=t*a[i][i];

printf("%lf\t",t);

}

printf("the determinant is:\n%lf",t);

}

void main()

{

input();

gauss();

determinant();

}

输出结果：

三、总结

程序编的太长如果把思路全部写出来，时间花的太多，因此具体总结在此就不写了。