九点差分格式.docx
《九点差分格式.docx》由会员分享，可在线阅读，更多相关《九点差分格式.docx（10页珍藏版）》请在冰豆网上搜索。
九点差分格式
二阶椭圆型方程的差分格式
题目：
九点差分格式
1.考虑问题
考虑Poisson方程:
（1）
G是xy平面上一有界区域，其边界
为分段光滑曲线。
在
上
满足下列边值条件之一：
（第一边值条件、强制边值条件），
（第二边值条件），
（第三边值条件、混合边值条件），
，
，
，
，及
都是连续函数，
。
2.九点差分格式
2.1公式推导
因为有
将（3）、（4）相加得：
又因为
因此
舍去截断误差，便得到逼近Poisson方程的九点差分格式：
其截断误差的阶为
化其形式为
=
同样将上述方程化成如下形式：
其中：
另外，
。
我们用
来表示数值解满足的线性方程组，则有：
（具有三对角矩阵的特征），其中：
，
，
，
，
。
2.2差分方程数值例子的求解
考虑问题
其精确解
2.3差分方程数值例子求解
考虑问题
其精确解
2.4求解图像
3.附录
f（x,y）函数
functionf=ff（x,y）
f=（pi*pi-1）*exp（x）*sin（pi*y）;
u（x,y）函数
functionf=fu（x,y）
f=sin（pi*y）*exp（x）;
clear;
clc;
formatshort;
%数据准备
M=40;%x轴划分M份
N=40;%y轴划分N份
x0=0;%x起点
xn=2;%x终点
y0=0;%y起点
yn=1;%y终点
Hx=（xn-x0）/（M+1）;
Hy=（yn-y0）/（N+1）;
%对自变量进行赋值
fori=1:
M
x（i）=x0+i*Hx;
ux0（i）=fu（x（i）,0）;
ux1（i）=fu（x（i）,1）;
end
forj=1:
N
y（j）=y0+j*Hy;
u0y（j）=fu（0,y（j））;
u2y（j）=fu（2,y（j））;
end
%a是u（i-1,j）系数b是u（i,j-1）系数c是u（i,j）系数d是u（i,j+1）系数e是u（i+1,j）系数
KK=-（Hx*Hx+Hy*Hy）/（12*Hx*Hx*Hy*Hy）;
a=KK;
b=-1/（Hx*Hx）-2*KK;
c=KK;
d=-1/（Hy*Hy）-2*KK;
e=2/（Hx*Hx）+2/（Hy*Hy）+4*KK;
f1=-1/（Hy*Hy）-2*KK;
g=KK;
h=-1/（Hx*Hx）-2*KK;
m=KK;
A=diag（ones（1,N）*e）+diag（ones（1,N-1）*f1,1）+diag（ones（1,N-1）*d,-1）;
I=diag（ones（1,N）*h）+diag（ones（1,N-1）*m,1）+diag（ones（1,N-1）*g,-1）;
J=diag（ones（1,N）*b）+diag（ones（1,N-1）*c,1）+diag（ones（1,N-1）*a,-1）;
maxA=blkdiag（kron（eye（M）,A））;
maxB=blkdiag（kron（diag（ones（1,M-1）,1）,I））;
maxC=blkdiag（kron（diag（ones（1,M-1）,-1）,J））;
max=maxA+maxB+maxC;
forii=1:
M
forjj=1:
N
f（（ii-1）*N+jj）=ff（x（ii）,y（jj））+1/12*（Hx*Hx*uxx（x（ii）,y（jj））+Hy*Hy*uyy（x（ii）,y（jj）））;
if（ii==1）
if（jj==1）
f（（ii-1）*N+jj）=f（（ii-1）*N+jj）-a*fu（x0,y0）-b*u0y（jj）-c*u0y（jj+1）-d*ux0（ii）-g*ux0（ii+1）;
end
if（jj==N）
f（（ii-1）*N+jj）=f（（ii-1）*N+jj）-a*u0y（jj-1）-b*u0y（jj）-c*fu（x0,yn）-f1*ux1（ii）-m*ux1（ii+1）;
%jj
end
if（jj~=1&jj~=N）
f（（ii-1）*N+jj）=f（（ii-1）*N+jj）-a*u0y（jj-1）-b*u0y（jj）-c*u0y（jj+1）;
end
end
if（ii==M）
if（jj==1）
f（（ii-1）*N+jj）=f（（ii-1）*N+jj）-a*ux0（ii-1）-d*ux0（ii）-g*fu（xn,y0）-h*u2y（jj）-m*u2y（jj+1）;
end
if（jj==N）
f（（ii-1）*N+jj）=f（（ii-1）*N+jj）-c*ux1（ii-1）-f1*ux1（ii）-g*u2y（jj-1）-h*u2y（jj）-m*fu（xn,yn）;
end
if（jj~=1&jj~=N）
f（（ii-1）*N+jj）=f（（ii-1）*N+jj）-g*u2y（jj-1）-h*u2y（jj）-m*u2y（jj+1）;
end
end
if（ii~=1&ii~=M）
if（jj==1）
f（（ii-1）*N+jj）=f（（ii-1）*N+jj）-a*ux0（ii-1）-d*ux0（ii）-g*ux0（ii+1）;
end
if（jj==N）
f（（ii-1）*N+jj）=f（（ii-1）*N+jj）-c*ux1（ii-1）-f1*ux1（ii）-m*ux1（ii+1）;
end
end
end
end
jg=（inv（max）*f'）';
jg;
jjgg=reshape（jg,M,N）;
foriii=1:
M
forjjj=1:
N
JJGG（iii,jjj）=fu（x（iii）,y（jjj））;
end
end
[X,Y]=meshgrid（x,y）;
subplot（1,2,1）;
surf（X,Y,jjgg）
title（'九点格式'）;
subplot（1,2,2）;
surf（X,Y,JJGG'）
title（'u（x,y）=sin（pi*y）*exp（x）'）;
升级会员 