九点差分格式.docx
- 文档编号:27560950
- 上传时间:2023-07-02
- 格式:DOCX
- 页数:10
- 大小:247.77KB
九点差分格式.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)');
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 九点差分 格式
![提示](https://static.bdocx.com/images/bang_tan.gif)