Possion方程的差分法求解.docx
- 文档编号:10402803
- 上传时间:2023-02-11
- 格式:DOCX
- 页数:12
- 大小:555.05KB
Possion方程的差分法求解.docx
《Possion方程的差分法求解.docx》由会员分享,可在线阅读,更多相关《Possion方程的差分法求解.docx(12页珍藏版)》请在冰豆网上搜索。
Possion方程的差分法求解
Possion方程的差分方法
课程名称:
微分方程数值解
学生姓名:
张弘
一、问题描述
二、问题分析
I.
偏微分方程的数值解法主要有有限差分法和Galerkin有限元法。
用差分法和有限元法将连续问题离散化的步骤是:
1、对求解区域做网格剖分用有限个网格节点代替连续区域。
2、微分算子离散化。
3、把微分方程的定解问题化为线性代数方程组的求解问题。
差分法和有限元方法的主要区别是离散化的第二步,差分法从定解问题的微分或积分形式出发,用数值微商或数值积分公式到处相应的线性代数方程组,有限元法从定解问题的变分形式出发,用Ritz-Galerkin法导出相应的线性代数方程组。
差分法的基本问题有:
(1)对求解区域做网格剖分
一维情形为把区间分为等距或不等距的区间,二维情形是把区域分割成均匀或不均匀的矩形,其边与坐标轴平行,也可分割成小三角形或凸四边形。
(2)构造逼近微分方程定解问题的差分格式
有两种构造差分格式的方法:
直接差分化法和有限体积法。
(3)差分解的存在唯一性,收敛性和稳定性研究
(4)差分方程的解法:
逐次超松弛法,交替方向法,共轭梯度法。
II
(1)由题目可知,本题需要考虑矩形网的五点差分格式。
题目给出的为第一边值条件,
将十字形图形的中心放于坐标原点处,如下图所示:
O
S1
G1
S2
由图形可知,所给区域可以看出是由8个梯形G1通过旋转、翻转拼接而成。
因此为了方便计算、减少计算量,只针对G1进行网格剖分,用5点差分格式进行求解。
但是由于G1是直角梯形,进行网格剖分时会出现x轴方向网格点个数不同的现象,不利于有差分形成的系数矩阵的生成,所以将三角形S1和梯形G1合在一起形成一个矩形,其区域为[0,3/2]×[0,1/2]。
如果采用等距差分,并且有hx=hy=h,设步长为h,
x=0:
h:
3/2;
y=0:
h:
1/2;
nx=length(x)-1;%为所求区域中x轴上内点的个数
ny=length(y)-1;%为所求区域中y轴上内点的个数
对于原来的梯形G1来说,有网格内点nx*ny-(ny^2-my)/2
对于矩形区域S1+G1而言,网格内点数为nx*ny,所以在后面的程序中要比在梯形G1中多计算了(ny^2-my)/2个点的函数值,但对程序效率的影响并不是很大,可以接受。
以下具体说明只需在G1上求解poisson方程的原因
所求方程为:
设直线l是经过原点o的任意一条直线,其方程为:
y=kx
设p(x,y)是区域内任一点,则其关于l对称的点为p'(s,t)
S=((k-1)^2+|2y)/(k^2+1)t=(2kx+(k^2-1)y)/(k^2+1)
则
同理可得:
将u(s,t)代替u(x,y)得:
Uxx+uyy=uss+utt=1
其依然满足上述方程。
令θ=arctan(k)点p的横坐标x=rcos(α)y=rsin(α)
则p关于直线l的对称点为p'(rcos(2θ-α),rsin(2θ-α)
由上述证明可知u(p)=u(p')。
由θ和p点的任意性知,对于函数u图像上的任意一点p,其关于任意一条经过原点直线l的对称点p'都在u的图像上,即u(α+δ)=u(α),即u关于原点是旋转对称的。
当l为x轴时,有u(x,y)=u(x,-y)
l为y轴时,u(x,y)=u(-x,y)
坐标轴旋转不改变方程的形式,所以函数在直角坐标系中u关于原点是旋转对称的,又所求十字形区域关于x,y轴是轴对称和关于原点中心对称的,因此可通过直求解区域G1,就可以知道函数在整个十字形区域的图像。
(2)对S1+G1形成的矩形进行正方形网格剖分,则求解Poisson方程的差分格式和化为如下形式:
对于正则内点其差分如下:
-Δuij=1/h^2*(-u(i,j+1)-u(i-1,j)+4u(i,j)-u(i+1,j)-u(i,j-1))=fij
(1)
对于矩形区域S1+G1,nx=(xb-xa)/h-1ny=(yb-ya)/h-1
按从左到右,从下到上的次序排列网点(ij):
j=1,1≤i≤nx;j=2,1≤i≤nx;,…;j=ny,1≤i≤nx,定义向量
Uh=[u11,u21…,unx-1;…;u1,nx-1,u2,nx-1,…,uny-1,nx-1]T
利用边界条件可以将
(1)式写成如下形式:
其中A=
为nx*ny阶矩阵,I为nx阶单位矩阵,B为nx阶单位矩阵。
B=
右端向量g的元素,依赖于边值a和右端项f,显然A是对称正定矩阵,也是稀疏矩阵因此可以采用逐次超松弛法,共轭梯度法和交替方向法莱求解,但为了方便程序设计采用了matlab的‘\'运算符来求解u。
针对本题的Poisson方程,对S1+G1形成的矩形网格的五点差分的具体分析。
对S1+G1形成的矩形五点差分的具体分析。
红色线条围成的区域为G1,L为红色斜边,S1为L左侧的三角形,S2为L右侧的三角形。
由对称性知,S1和S2中关于L相互对称的网格点其上U的函数值是相同的。
按从左到右,从下到上的次序排列网点(ij)。
(1)
对于三角形S1中的网点有u(i,j),ny≥i>j≥1,有u(i,j)-u(j,i)=0
所以S1中网点的差分就为:
u(i,j)-u(j,i)=0
(2)由于原点o点的特殊性,其邻点中u(1,2)=u(1,-1)=u(-1,1)=u(2,1)
所以其差分为:
u(1,1)-4u(1,2)=h^2*f(i,j)
程序中语句:
A(1:
nx,nx+1:
2*nx)=2*I;和A(1,2)=-2;
保障了上面差分方程的系数。
(3)对于网格上最下面一行上除了原点o的所有正则网格内点,由对称性得u(1,i)的邻点中
u(1,i-1)=u(1,i+1)
所以其差分为:
4u(i,j)-u(i-1,j)-u(i+1,j)-2*u(i,j+1)=h^2*f(i,j)
对于右下角的非正则内点u(1,nx)
其差分为:
4u(i,j)-u(i-1,j)-2*u(i,j+1)=h^2*f(i,j)
程序中的相关语句为:
A(1:
nx,nx+1:
2*nx)=2*I;
A(nx+1:
2*nx,nx+1:
2*nx)=B;
(4)对于G1中的所有正则内点,其差分为:
4u(i,j)-u(i-1,j)-u(i+1,j)-u(i,j-1)-u(i,j+1)=h^2*f(i,j)
程序中相关语句:
A(i*nx+1:
(i+1)*nx,i*nx+1:
(i+1)*nx)=B;
A((i-1)*nx+1:
i*nx,i*nx+1:
(i+1)*nx)=I;
A(i*nx+1:
(i+1)*nx,(i-1)*nx+1:
i*nx)=I;
(5)对于网格最后一列的所有非正则内点u(:
nx),其差分为:
4u(nx,j)-u(nx,j-1)-u(nx,j+1)-u(nx-1,j)=h^2*f(i,j)
(6)在求出了所有的内点后,对于剩下的边界点赋值有:
U(ny+1,ny+1:
nx)=0%最上面一行上的边界点
U(1:
ny+1,nx+1)=0%最右侧一列的边界点
u(ny+1,1:
ny)=u(1:
ny,ny+1);%三角形S1和S2边界上的网点,它们是S1的边界点,但是整个求解区域的内点。
三、程序设计及分析
functionpoisson5spot(h)
ifnargin<1%默认步长h=
h=;
end
x=0:
h:
3/2;
y=0:
h:
1/2;
nx=length(x)-1;%取区域的内点
ny=length(y)-1;
%===========构造矩阵B==========
B=eye(nx,nx)*4;
fori=1:
nx-1
B(i,i+1)=-1;
end
fori=2:
nx
B(i,i-1)=-1;
end
I=-eye(nx,nx);
%===========构造系数矩阵A========
A=zeros(nx*ny,nx*ny);
A(1:
nx,1:
nx)=B;
%由区域的对称性,正方形网格最下面一行的差分形式
%改为4u(i,j)-u(i-1,j)-u(i+1,j)-2*u(i,j+1)
%因为u(i,j+1)=u(i,j-1)
A(1:
nx,nx+1:
2*nx)=2*I;
A(nx+1:
2*nx,nx+1:
2*nx)=B;
%矩形网格左下角第一个点的差分形式改为:
%4u(i,j)-u(i-1,j)-u(i+1,j)-u(i,j+1)-u(i,j-1)
%=4u(i,j)-2u(i,j+1)-2u(i+1,j)
A(1,2)=-2;
%为了方便,本程序往梯形中增加了一个三角形,以方便编程
A(nx+1:
2*nx,1:
nx)=I;
fori=2:
ny-1
A(i*nx+1:
(i+1)*nx,i*nx+1:
(i+1)*nx)=B;
A((i-1)*nx+1:
i*nx,i*nx+1:
(i+1)*nx)=I;
A(i*nx+1:
(i+1)*nx,(i-1)*nx+1:
i*nx)=I;
end
%==========================
%bi=h*h*f;右端向量
b=h*h*ones(nx*ny,1);
%由于对称性,左侧三角形中有u(i,j)-u(j,i)=0i>j
%因此令A(i,j)=1A(j,i)=-1
%所以在本程序中多计算了(ny^2-my)/2个点的函数值
%但对程序的影响并不是很大
fori=2:
ny
A((i-1)*nx+1:
(i-1)*nx+i-1,:
)=0;
forj=1:
i-1
A((i-1)*nx+j,(i-1)*nx+j)=1;
A((i-1)*nx+j,(j-1)*nx+i)=-1;
b((i-1)*nx+j)=0;
end
end
x=A\b;%为了方便采用左除求解网格点上的函数值
%x=gmres(A,b,100);
%按顺序将x赋值给u
u=zeros(ny+1,nx+1);
fori=1:
ny
forj=1:
nx
u(i,j)=x((i-1)*nx+j);
end
end
%根据对称性,给网格最上面一行的点赋值
u(ny+1,1:
ny)=u(1:
ny,ny+1);
%=============作Poisson方程在区域上的图形===========
[x,y]=meshgrid(0:
h:
3/2,0:
h:
1/2);
holdon
surf(x,y,u)%11第一象限的第一块区域,下面的以此类推
surf(y,x,u)%12
surf(-y,x,u)%21
surf(-x,y,u);%22
surf(-x,-y,u)%31
surf(-y,-x,u)%32
surf(y,-x,u)%41
surf(x,-y,u);%42
四、实验结果
1.在区域G1上求解后的图形显示:
图形表示了在S1+G1区域上的函数图像,而不是单纯的G1上的函数图像。
2通过拼接后图形显示如下:
由上图可知,除了边界点外网格点上的函数值都有u(i,j)>0.
Lhuij>0,对任意(xi,yj)∈Gh,uij不可能在内点取得负的极小,与极值定理相符合。
3、翻转后图形如下:
4网格间距h=时得到的图形:
五、实验分析
本次实验将求解区域先利用对称性将其缩小为原区域的1/8,又为了方便系数矩阵的选取给G1添加了三角形S1,减少了运算量。
在实验中通过更改h,发现当h=时,使用matlab的‘\‘符号来求解,会发生内存溢出的现象,因此此时系数矩阵A为500*1500=750000,但是如果采用稳定双共轭梯度法bicgstab和预条件的再开始gmres算法则可以继续求解。
6、参考文献
[1]李荣华,刘播.微分方程数值方法.
[2]胡建伟,汤怀民,微分方程的数值方法.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Possion 方程 差分法 求解