![基于高斯消元法的三对角矩阵LU分解_第1页](http://file4.renrendoc.com/view/5e6516e6443d7130f06c8cf1e03039db/5e6516e6443d7130f06c8cf1e03039db1.gif)
![基于高斯消元法的三对角矩阵LU分解_第2页](http://file4.renrendoc.com/view/5e6516e6443d7130f06c8cf1e03039db/5e6516e6443d7130f06c8cf1e03039db2.gif)
![基于高斯消元法的三对角矩阵LU分解_第3页](http://file4.renrendoc.com/view/5e6516e6443d7130f06c8cf1e03039db/5e6516e6443d7130f06c8cf1e03039db3.gif)
![基于高斯消元法的三对角矩阵LU分解_第4页](http://file4.renrendoc.com/view/5e6516e6443d7130f06c8cf1e03039db/5e6516e6443d7130f06c8cf1e03039db4.gif)
![基于高斯消元法的三对角矩阵LU分解_第5页](http://file4.renrendoc.com/view/5e6516e6443d7130f06c8cf1e03039db/5e6516e6443d7130f06c8cf1e03039db5.gif)
版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
本文格式为Word版,下载可任意编辑——基于高斯消元法的三对角矩阵LU分解三对角与块三对角方程组课程设计
一、基于高斯消元法的三对角方程组求解
三对角矩阵是一类重要的特别矩阵,在数学计算和工程计算中有广泛应用。例如,二阶常微分方程边值问题数值求解,一维热传导方程数值求解,以及三次样条函数计算等都会涉及到三对角方程组求解。由于三对角矩阵的稀疏性质,用直接法求解三对角方程组的算法效率较高,很有实用价值。
考虑n阶三对角矩阵和n维向量
?f1???1?1??f??????122?,f??2?A=???????????????n?1n???fn?求解方程组Ax=f的高斯消元法的程序如下
functionf=triGauss(gama,alpha,bata,f)
%SolvingTriDiag(gama,alpha,bata)systemsbyGaussmethodn=length(alpha);fork=1:n-1
m=gama(k)/alpha(k);
alpha(k+1)=alpha(k+1)-m*bata(k);f(k+1)=f(k+1)-m*f(k);end
f(n)=f(n)/alpha(n);fork=n-1:-1:1
f(k)=(f(k)-bata(k)*f(k+1))/alpha(k);end
由程序知,对于n阶三对角方程组,高斯消元法只用到5n–4次乘法和除法。例1.求二阶常微分方程边值问题
??y???y?x,x?(0,1)?y(0)?0,y(1)?0?sinhx数值解,并与解析解:y(x)?x?作对比。
sinh1解:对正整数n,取h=1/(n+1),令xj=jh,(j=0,1,2,……,n+1),yj≈y(xj)。由
yj?1?2yj?yj?1y?j??2h得差分方程
?整理,得
yj?1?2yj?yj?1h2?yj?xj
–yj-1+(2+h2)yj–yj+1=h2xj(j=1,…,n)
根据边界条件,有
y0=0,yn+1=0
形成三对角方程组
?2?h2???1??????12?h2??1???12?h2?1?x1???y1?????y?x?2???2?2??h???????????x?1??yn?1??n?1????x?2?h2???yn??n?取n=9,得数值解和解析解的最大误差为:4.4146e-005,将数值解和解析解数据绘图如下
图1.数值解与解析解比较
程序如下
n=input('inputn:=');h=1/(n+1);x=0:h:1;
ux=x-sinh(x)./sinh(1);alpha=(2+h*h)*ones(n,1);bata=-ones(n-1,1);gama=bata;f(1:n)=h*h*x(2:n+1);
uk=triGauss(gama,alpha,bata,f);Uk=[0,uk,0];
error=max(abs(ux-Uk))plot(x,ux,x,Uk,'ro')
取不同的n,运行程序分别计算试验数据如下9192939n4.4146e-0051.1047e-0054.9106e-0062.7624e-006error数据说明,随着结点增加,数值解的误差变小。
二、迭代法求解三对角方程组
数值分析介绍的三种迭代格式用于例1中三对角方程组:1.JACOBI迭代
u(jk?1)?u(jk?1)?12(k)(k)[hx?u?ujj?1j?1]22?h1?1)(k)[h2xj?u(jk?1?uj?1]22?h2.SEIDEL迭代
3.SOR迭代
u(jk?1)?(1??)u(jk)??2?h2(k?1)(k)[hx?u?ujj?1j?1]2由于JACOBI迭代矩阵的谱半径
?(BJ)?SEIDEL迭代矩阵的谱半径
2cos?h
2?h2222)cos?h22?h?(BG?S)?(根据,的结果,取SOR迭代松驰因子为
??SOR迭代矩阵的谱半径
21?1??(BJ)2
?(BSOR)???1?1?1??(BJ)21?1??(BJ)2
试验数据如下(迭代精度要求10-8)nJACOBISEIDELSOR迭代数误差迭代数误差迭代次数误差5880.1191e-3470.1190e-3180.1190e-392300.4432e-41230.4422e-4290.4415e-4198290.1173e-44420.1137e-4560.1106e-43929130.5635e-515610.4176e-51070.2789e-5数据说明,三种迭代法计算出的数值解误差一致。但对三对角方程组而言,尽管迭代法程序比较简单,但迭代法的效率不如直接法。程序如下
n=input('inputn:=');h=1/(n+1);x=0:h:1;y=x-sinh(x)/sinh(1);hh=h*h;rr=2+hh;u0=zeros(size(x));u=u0;k1=0;error=1.0;%Jacobi迭代whileerror>1.0e-8error=0;forj=2:n+1
temp=(hh*x(j)+u(j-1)+u(j+1))/rr;error=max(error,abs(temp-u(j)));v(j)=temp;end
u=[v,0];k1=k1+1;end
error1=max(abs(u-y));k1
u=u0;error=1;k2=0;%Seidel迭代whileerror>1.0e-8error=0;forj=2:n+1temp=u(j);
u(j)=(hh*x(j)+u(j-1)+u(j+1))/rr;error=max(abs(temp-u(j)),error);end
k2=k2+1;end
error2=max(abs(u-y));k2
u=u0;error=1;k3=0;%SOR迭代ro=2*cos(h*pi)/rr;
omiga=2/(1+sqrt(1-ro*ro));omiga1=1-omiga;whileerror>1.0e-8error=0;forj=2:n+1temp=u(j);
u(j)=omiga1*temp+omiga*(hh*x(j)+u(j-1)+u(j+1))/rr;error=max(abs(temp-u(j)),error);end
k3=k3+1;end
error3=max(abs(u-y));k3
[error1,error2,error3][ro,ro*ro,-omiga1]
三、块三对角方程组的迭代方法
拉普拉斯(Laplace)是法国著名数学家和天文学家。他用数学方法证明白行星的轨道大小只有周期性变化,这就是著名拉普拉斯的定理。1796年,他发表《宇宙体系论》,因研究太阳系稳定性的动力学问题被誉为法国的牛顿和天体力学之父。拉普拉斯在数学和物理学方面也有重要贡献,以他的名字命名的拉普拉斯变换和拉普拉斯方程,在科学技术的各个领域有着广泛的应用。对于二维拉普拉斯方程问题,方程的解是二元函数。对给定的边界条件,求数值解需要解块三对角方程组。
例2.正方形区域上狄利克莱边值问题
?uxx?uyy?0,0?x,y?1??u(0,y)?u(x,0)?u(x,1)?0?u(1,y)?sin?y?sh?xsin?y。sh?取正整数n,记h=1/(n+1),对区域做网格剖分:
xi=ih(i=0,1,……,n+1)yj=jh(j=0,1,……,n+1)
在点(xi,yj)处记uij=u(xi,yj)。离散拉普拉斯方程,可得常见的的五点处差分格式
ui-1,j+ui,j-1–4uij+ui+1,j+ui,j+1=0(i,j=1,…,n)
u0,j=0,un,j=sin(πj/n)(j=1,…,n)ui,0=0,ui,n=0(i=1,…,n)
确凿解:u(x,y)?该方程组的系数矩阵是块三对角矩阵,与三对角矩阵相比,依旧具有稀疏性。但问题的规模比较大了。对于n=4,其矩阵是16阶方程,矩阵非零元分布如下图所示
图2.块三对角矩阵非零元分布图
1.五点差分格式的Seidel迭代方法
1k?1)(k?1)(k)(k)ui(,kj?1)?[ui(?1,j?ui,j?1?ui?1,j?ui,j?1],(i,j=0,1,…,N)
42.SOR迭代算法:
(k)ui(,kj?1)?(1??)uij??4k?1)(k?1)(k)(k)[ui(?1,j?ui,j?1?ui?1,j?ui,j?1],(i,j=0,1,…,N)
其中,最正确松驰因子为
??2
1?sin?h6024291260.48407.3660e-005五点差分格式seidle迭代试验数据表(误差限要求:1e-008)2102202402结点数n1826062077迭代次数0.29704.1158.531CPU时间0.00236.4274e-0041.6814e-004误差超松驰迭代法,取最正确松驰因子:??2
1?sin?h五点差分格式SOR迭代试验数据表(误差限要求:1e-008)2102202402602结点数n4074139201迭代次数0.110.65604.953015.5940CPU时间0.00236.4306e-0041.6944e-0047.6600e-005误差对比两种迭代方法,知SOR迭代无论从迭代次数还是从计算机运行时间比较都优于seidle迭代法。但两种方法的数值解误差几乎是一样的。迭代计算绘图显示如下
图3.二维拉普拉斯方程边值问题解图形
由于边界条件中,右侧边界条件非零,而其余边界条件均为零,故有此图形。图中红色区域表
示对应的函数值数值大,蓝色区域部分则说明对应的函数值数值小。
附:MATLAB程序
1.五点差分格式的seidel迭代法程序:
n=input('n=');
h=1/(n+1);x=0:h:1;y=x';u=zeros(n+2,n+2);u(:,n+2)=sin(pi*y);k=0;er=1;t=cputime;
whileer>1e-008er=0;k=k+1;fori=2:n+1forj=2:n+1
s=(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))/4;er=max(er,abs(s-u(i,j)));u(i,j)=s;endendendk
times=cputime-tpcolor(u)
w=sin(pi*y)*sinh(pi*x)/sinh(pi);error=max(max(abs(w-u)))
2.五点差分格式SOR迭代程序
n=input('n=');
h=1/(n+1);x=0:h:1;y=x';u=zeros(n+2,n+2);u(:,n+2)=sin(pi*y);w=2/(1+sin(pi*h));w1=1-w;k=0;er=1;t=cputime;
whileer>1e-008er=0;k=k+1;fori=2:n+1forj=2:n+1s=u(i,j);
u(i,j)=w1*s+w*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))/4;er=max(er,abs(s-u(i,j)));
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 商铺门头分租合同协议书
- 个人借款分期合同协议书
- 公司签不平等合同协议书
- 出售防水器材合同协议书
- 合同协议书购销合同模板
- 合同协议书模板码
- 主播网签合同协议书
- 园林苗木订购合同协议书
- MMG-23600-TPU水凝胶屏幕保护膜市场调研报告全球行业规模展望2024-2030 Sample
- 酒店实美学试题库(二)
- 2024年秋季新学期开学国旗下讲话稿(二篇)
- 大规模设备更新 消费品以旧换新 政策汇编(国家)
- 酒店入股合同协议书2024年
- 2024-2030年连锁超市行业市场发展现状及并购重组策略与投融资研究报告
- DB32T 4830-2024蛋鸡养殖场健康生产规范
- 新质农业劳动力培育对农作物生产的影响
- 2024巴黎奥运会秋季开学第一课主题班会
- GB/T 20878-2024不锈钢牌号及化学成分
- 2024年河南客车驾驶员从业资格证考试题库
- 外研版(三起)(2024)三年级上册英语Unit 1 Let's be friends 教学设计
- 《财务大数据分析》 课件 项目1 财务大数据分析认知
评论
0/150
提交评论