5点差分格式的Matlab程序.docx
- 文档编号:7297017
- 上传时间:2023-01-22
- 格式:DOCX
- 页数:12
- 大小:35.17KB
5点差分格式的Matlab程序.docx
《5点差分格式的Matlab程序.docx》由会员分享,可在线阅读,更多相关《5点差分格式的Matlab程序.docx(12页珍藏版)》请在冰豆网上搜索。
5点差分格式的Matlab程序
5点差分格式的Matlab程序
三对角与块三对角方程组课程设计
一、基于高斯消元法的三对角方程组求解
三对角矩阵是一类重要的特殊矩阵,在数学计算和工程计算中有广泛应用。
例如,二阶常微分方程边值问题数值求解,一维热传导方程数值求解,以及三次样条函数计算等都会涉及到三对角方程组求解。
由于三对角矩阵的稀疏性质,用直接法求解三对角方程组的算法效率较高,很有实用价值。
考虑n阶三对角矩阵和n维向量
,f,,,,111,,,,,,,f1222,,,,A=,,f,,,,?
?
?
?
,,,,,fn,1nn,,,,
求解方程组Ax=f的高斯消元法的程序如下
functionf=triGauss(gama,alpha,bata,f)
%SolvingTriDiag(gama,alpha,bata)systemsbyGaussmethod
n=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,
sinhxy(x),x,数值解,并与解析解:
作对比。
sinh1
解:
对正整数n,取h=1/(n+1),令x=jh,(j=0,1,2,……,n+1),y?
y(x)。
由jjj
y,2y,y11j,jj,,,y,j2h
得差分方程
y,y,y2jjj,1,1,,y,xjj2h
整理,得22–y+(2+h)y–y=hx(j=1,…,n)j-1jj+1j
根据边界条件,有
y=0,y=00n+1
形成三对角方程组
2yx,,,,,,2,h,111,,,,,,2yx,12,h,122,,,,,,2,,,,,,?
h?
?
?
?
,,,,,2yx,12,h,1n,1n,1,,,,,,2,,,,,,yx12h,,nn,,,,,,取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,运行程序分别计算实验数据如下
9192939n
4.4146e-0051.1047e-0054.9106e-0062.7624e-006error
数据表明,随着结点增加,数值解的误差变小。
二、迭代法求解三对角方程组
数值分析介绍的三种迭代格式用于例1中三对角方程组:
1(JACOBI迭代
1,(k1)2(k)(k)u,[hx,u,u],,jjj1j122,h
2(SEIDEL迭代
1,,(k1)2(k1)(k)u,[hx,u,u],,jjj1j122,h
3(SOR迭代
,,(k1)(k)2(k1)(k)u,(1,)u,[hx,u,u],,,jjjj1j122,h
由于JACOBI迭代矩阵的谱半径
2,,(B)cos,hJ2,2h
SEIDEL迭代矩阵的谱半径
222,(B),()cos,hG,S22,h
根据,的结果,取SOR迭代松驰因子为
2,,21,1,,(B)J
SOR迭代矩阵的谱半径
2,1,1,(B)J,,(B),,1,SOR21,1,,(B)J-8实验数据如下(迭代精度要求10)
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-8
error=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-8
error=0;
forj=2:
n+1
temp=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-8
error=0;
forj=2:
n+1
temp=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(正方形区域上狄利克莱边值问题
0,0,1u,u,,xy,,xxyy,(0,)(,0)(,1)0uy,ux,ux,,
(1,)sin,uy,y,
shx准确解:
。
u(x,y),sin,ysh,
取正整数n,记h=1/(n+1),对区域做网格剖分:
x=ih(i=0,1,……,n+1)i
y=jh(j=0,1,……,n+1)j
在点(x,y)处记u=u(x,y)。
离散拉普拉斯方程,可得常见的的五点处差分格式ijijij
u+u–4u+u+u=0i-1,ji,j-1iji+1,ji,j+1
(i,j=1,…,n)
u=0,u=sin(πj/n)(j=1,…,n)0,jn,j
u=0,u=0(i=1,…,n)i,0i,n
该方程组的系数矩阵是块三对角矩阵,与三对角矩阵相比,仍然具有稀疏性。
但问题的规模比较大了。
对于n=4,其矩阵是16阶方程,矩阵非零元分布如下图所示
图2.块三对角矩阵非零元分布图1(五点差分格式的Seidel迭代方法
1,,,(k1)(k1)(k1)(k)(k),(i,j=0,1,…,N)u,[u,u,u,u],,,,i,ji1,ji,j1i1,ji,j14
2(SOR迭代算法:
,,,(k1)(k)(k1)(k1)(k)(k),(i,j=0,1,…,N)u,(1,)u,[u,u,u,u],,,,,i,jiji1,ji,j1i1,ji,j14
其中,最佳松驰因子为
2,,1,sin,h五点差分格式seidle迭代实验数据表(误差限要求:
1e-008)2222210204060结点数n
18260620774291迭代次数
0.29704.1158.531260.4840CPU时间
0.00236.4274e-0041.6814e-0047.3660e-005误差
2,超松驰迭代法,取最佳松驰因子:
1,sin,h
五点差分格式SOR迭代实验数据表(误差限要求:
1e-008)2222210204060结点数n
4074139201迭代次数
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-008
er=0;k=k+1;
fori=2:
n+1
forj=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;
end
end
end
k
times=cputime-t
pcolor(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-008
er=0;k=k+1;
fori=2:
n+1
forj=2:
n+1
s=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)));
end
end
end
k
times=cputime-t
%surf(u)
%figure,pcolor(u)
z=sin(pi*y)*sinh(pi*x)/sinh(pi);
error=max(max(abs(z-u)))
实验题
1(求下列二阶常微分方程数值解
x,,,u(x),u(x),e(sinx,2cosx),0,x,,
u(0),0,u(,),0
x该问题解析解为:
u(x)=esinx
2(求正方形区域上的泊松方程数值解
222,0 边界条件 12u(x,0)=0,u(x,1),x,0 1,2,u(1,y),esiny,y,0 1,x2精确解: u(x,y),esin,y,(xy)2 利用五点有限差分法求解含Dirichlet条件的Poisson问题的近似解 function[u,x,y,error]=poissonfd(a,b,c,d,nx,ny,fun,bound,uex,varargin)%POISSONFDtwo-dimensionalPoissonsolver%[u,x,y]=POISSONFD(A,B,C,D,NX,NY,FUN,BOUND)solvesby%the5-pointsfinitedifferenceschemetheproblem%-LAQPL(U)=Funintehrectangle(A,C)X(B,D)with%DirichletboundaryconditionsU(X,Y)=BOUND(x,y)forany%(x,y)ontheboundaryoftherectangle.% %[U,X,Y,ERROR]=POISSONFD(A,C,B,D,NX,NY,FUN,BOUND,UEX)compu'%alsothemaximumnodalerrorERRORwithrespecttotheexact %solutionUEX. %FUN,BOUNDandUEXcanbeonlinefunctions.ifnargin==8 uex=inline('0','x','y'); end nx=nx+1; ny=ny+1; hx=(b-a)/nx; hy=(d-c)/ny; nx1=nx+1;hx2=hx^2; hy2=hy^2; kii=2/hx2+2/hy2;kix=-1/hx2;kiy=-1/hy2;dim=(nx+1)*(ny+1); k=speye(dim,dim); rhs=zeros(dim,I); y=c; form=2: ny x=a;y=y+hy; forn=2: nx i=n+(m-1)*(nx+1); x=x+hx; rhs(i)=feval_r(fun,x,y,varargin{: });K(i,i)=kii; k(i,i-1)=kix; k(i,i+1)=kix; k(i,i+nx1)=kiy; K(i,i-nx1)=kiy; end end rhsl=zeros(dim,i); x=[a: hx: b]; rhsl(1: nxl)=feval_r(bound,x,c,varargin{: }); rhsl(dim-nx: dim)=feval_r(bound,x,d,varargin{;}); y=[c: hy: d]; rhsl(1: nx1: dim-nx)=feval_r(bound,a,y,varargin{: }); rhsl(nx1: nx1: dim)=feval_r(bound,b,y,varargin{: }); rhs=rhs=K*rhsl; nbound=[[1: nx1],[dim-nx: dim],[1: nx1: dim-nx],[nx1: nx1: dim]]; ninternal=setdiff([1: dim],nbound);K=K(ninternal,ninternal); Rhs=rhs(ninternal); Utemp=K\rhs; uh=rhsl; uh(ninternal)=utemp; k=1;y=c; forj=1: ny+1 x=a; fori=1: nx1 u(i,j)=uh(k); k=k+1; ue(i,j)=uh(k); k=k+1; ue(i,j)=feval_r(uex,x,y,varargin{: });x=x+hx; end y=y+hy; end x=[a: hx: b]; y=[c: hy: d]; ifnargout==4&nargin==8 warning('Exactsolutionnotavailable');error=[]l else error=max(max(abs(u-ue)))/max(max(abs(ue))); end end return feval函数如何理解 feval函数的最通常的应用是以下形式: feval(‘functionname’,parameter),举个简单的例子: 比如要计算sin (2),当然可以直接用命令y=sin (2);利用feval,还可以这 样来做: y,feval(‘sin’,2);另外这里的函数名字还可以是一个函数句柄,即 h=@sin; y=feval(h,2);或者直接写成y=feval(@sin,2); [y1,..,yn]=Feval_r(F,x1,...,xn)F是需要使用函数的函数名,或者句柄;xi是函数的参数,yi是函数的 返回值 举例: 假设需要调用的函数foo定义如下: functionx=foo(a,b) x=a*b; 若在main函数中用feval调用foo,可以有以下几种方式 1.result=feval_r('foo',3,15); 2.result=feval_r(@foo,3,16);%这里@foo即句柄 3.若调用的函数要作为main的参数,则 functionresult=main(f) result=feval_r(f,3,10); 然后调用main时将'foo'传入即可 >>main('foo');
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 点差分 格式 Matlab 程序
![提示](https://static.bdocx.com/images/bang_tan.gif)