crout分解法1.docx
- 文档编号:4155029
- 上传时间:2022-11-28
- 格式:DOCX
- 页数:9
- 大小:111.51KB
crout分解法1.docx
《crout分解法1.docx》由会员分享,可在线阅读,更多相关《crout分解法1.docx(9页珍藏版)》请在冰豆网上搜索。
crout分解法1
crout分解法
(1)
Crout分解法解线性方程组的算法及程序设计
【摘要】在自然科学和工程技术中很多问题的解决常常归结为解线性方程组,而方程组的系数矩阵大致分为两种,一种是低阶稠密矩阵(阶数不超过150),另一种是大型稀疏矩阵(矩阵阶数高且零元素较多).解低阶稠密矩阵和某些特殊形式的大型稀疏矩阵宜应用直接法,例如:
高斯消元法、矩阵三角分解法。
对于一般形式的系数矩阵,Crout分解法是一种有效的方法,本文就Grout分解法给出其计算公式的推导过程及算法程序。
【正文】
一、Crout方法解线性方程组的算法
给定线性方程组Ax=b,其中系数矩阵A=(aij)n×n为可逆的,x=(x1,x2,…,xn)T,b=(b1,b2,…bn)T为常数项列向量。
Crout分解后A=LU,L和U的结构为
(1)公式推导:
1由A=LU及矩阵乘积和相等概念,有
ai1=(li1,li2···,lii,0,···0)·
=li1(i=1,2,···n)
a1j=(l11,0,0,0···0)·
=l11·u1j(j=1,2,···,n)
得:
li1=ai1;u1j=a1j/l11(i,j=1,2,···,n)
2当i i=1,2,…n) ②解两个三角型方程组 由A=LU及Ax=b可得L(Ux)=b,令Y=(y1,y2,…,yn)T=Ux,则LY=b,于是求解Ax=b就被化为求解下三角型方程组LY=b及单位上三角型方程组Ux=Y。 a)先解下三角型方程LY=b。 由l11y1=b1,l21y1+l22y2=b2l11,……ln1y1+ln2y2+…+lnnyn=bn 所以y1=b1/l11,yi=(bi- )/lii,i=2,3,…,n b)再解单位上三角型方程组Ux=Y。 由UX=Y得x1+u12x2+…+u1nxn=y1,x2+…+u2nxn=y2 ………xn-1+un-1nxn=yn-1,xn=yn,利用回代解法可得方程组AX=b的解为 xn=yn,xi=yi- i=n-1,…,2,1 二、Crout方法解线性方程组的程序 (1)程序代码: #include"stdio.h" #include"math.h"//头文件 #defineN20//自定义N=20 intmain()//主函数 { inti,j,k; intsize; floata[N][N],l[N][N],u[N][N]; floatb[N],x[N],y[N];//定义变量 printf("\t\t\tCrout分解法解方程组\n"); printf("请输入方阵A的n: "); scanf("%d",&size); printf("\n"); printf("请输入方程组的系数: \n"); for(i=0;i { for(j=0;j { scanf("%f",&a[i][j]);//输入方程组系数矩阵a[][] } } printf("\n请输入方程组的y: \n"); for(i=0;i scanf("%f",&b[i]);//输入结果矩阵b[] printf("\n方阵A[][]为: \n"); for(i=0;i { for(j=0;j { printf("%f",a[i][j]);//输出a[][] } printf("\n"); } printf("\n方程组y为: \n"); for(i=0;i printf("%f",b[i]);//输出b[] printf("\n"); for(i=0;i { u[i][i]=1;//定初始值令u[i][i]=1 } for(i=0;i for(j=i+1;j { l[i][j]=0;//定初始值令l[i][j]=0 } for(j=0;j for(i=j+1;i { u[i][j]=0;//定初始值令u[i][j]=0 } l[0][0]=a[0][0]; for(i=1;i { l[i][0]=a[i][0];//计算第一行的l[][] u[0][i]=a[0][i]/l[0][0];//计算第一列的u[][] } for(i=1;i { for(j=1;j<=i;j++)//计算第2行到第size-1行的l[][] { l[i][j]=a[i][j]; for(k=0;k { l[i][j]=l[i][j]-l[i][k]*u[k][j]; } } printf("\n"); for(j=i+1;j { u[i][j]=a[i][j]; for(k=0;k<=i-1;k++) { u[i][j]=u[i][j]-l[i][k]*u[k][j]; } u[i][j]=u[i][j]/l[i][i]; } printf("\n"); } for(j=1;j { l[size-1][j]=a[size-1][j]; for(k=0;k<=j-1;k++) { l[size-1][j]=l[size-1][j]-l[size-1][k]*u[k][j]; } } printf("\n"); printf("输出矩阵L[i][j]\n"); for(i=0;i for(j=0;j { printf("%f",l[i][j]);printf("");//输出下三角矩阵l[][] } printf("\n"); } printf("输出矩阵U[i][j]\n"); for(i=0;i for(j=0;j { printf("%f",u[i][j]);printf("");//输出单位上三角矩阵u[][] } printf("\n"); } y[0]=b[0]/l[0][0];//给y[0]初始值 for(i=1;i { y[i]=b[i]; for(k=0;k<=i-1;k++) { y[i]=y[i]-l[i][k]*y[k];//计算公式 } y[i]=y[i]/l[i][i]; } printf("\n"); printf("y值: \n"); for(i=0;i printf("y[%d]=%f",i+1,y[i]);//输出y[i]的结果 printf("\n\n"); printf("x的值: \n"); x[size-1]=y[size-1];//给x[size-1]赋值 for(i=size-2;i>=0;i--) { x[i]=y[i]; for(k=i+1;k { x[i]=x[i]-u[i][k]*x[k];//计算x[] } } for(i=0;i { printf("x[%d]=%f\n",i+1,x[i]);//输出x[i]的结果 } } (2)程序应用
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- crout 解法