crout分解法1.docx
- 文档编号:15720704
- 上传时间:2023-07-07
- 格式:DOCX
- 页数:11
- 大小:113KB
crout分解法1.docx
《crout分解法1.docx》由会员分享,可在线阅读,更多相关《crout分解法1.docx(11页珍藏版)》请在冰点文库上搜索。
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的结构为
aij=(l11,li2,···lii,0···0)·
=
当i
j时
aij=(l11,li2,···lii,0···0)·
=
得:
lij=aij-
(j=2,3,···,i,i=2,3,···,n)
uij=(aij-
)/lii(i=2,3,···,n,j=i+1,i+2,···,n)
(2)具体算法如下:
①对A作LU分解:
由A=LU及矩阵的乘法原理可得:
lij=aij-
(j=1,2,…,i,i=1,2,…n)
uij=(aij-
)/lii(j=i+1,i+2,…,n,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 解法