偏微分方程上机实验报告.doc
- 格式:doc
- 大小:72.50 KB
- 文档页数:6
上机实验2:五点差分格式法
偏微分方程(Matlab )实验报告
——五点差分格式法 一、 实验题目
设G 是形如下图的十字形域,由五个相等的单位正方形组成,用五点差分格式求下列边值问题的数值解:
22
2
21,u u
G x y ∂∂+=-∂∂∂于u=0,于G
二、 实验原理
取定沿X 轴和Y 轴方向的步长1h 和2h ,()
12
22
1
2
h h h =+,作两族与坐
标轴平行的直线:x=i 1h ,y=j 2h ,,0,1,2,i j =±±
若(,i j x y )为正则内点,沿x,y 方向分别用二阶中心差商代替
xx yy u u 和则得
1,1,,1,1
2
212
22[
]i j ij i j
i j ij i j ij u u u u u u f h h +-+--+-+-+
=
特别取正方形网格:12h h h ==,则原差分方程可简化为
2
1,,11,,11()44
ij i j i j i j i j ij h u u u u u f --++-+++=
三、 实验程序
1)function uxy = EllIni2Uxl(x,y)
format long ;
uxy = 0;
2)function uxy = EllIni2Uxr(x,y)
format long;
uxy = y*(2-y);
3)function uxy = EllIni2Uyl(x,y)
format long;
uxy = 0;
4)function uxy = EllIni2Uyr(x,y)
format long;
if x < 1
uxy = x;
else
uxy = 2 - x;
end
5)function u = peEllip5(nx,minx,maxx,ny,miny,maxy) format long;
hx = (maxx-minx)/(nx-1);
hy = (maxy-miny)/(ny-1);
u0 = zeros(nx,ny);
for j=1:ny
u0(j,1) = EllIni2Uxl(minx,miny+(j-1)*hy);
u0(j,nx) = EllIni2Uxr(maxx,miny+(j-1)*hy);
end
for j=1:nx
u0(1,j) = EllIni2Uyl(minx+(j-1)*hx,miny);
u0(ny,j) = EllIni2Uyr(minx+(j-1)*hx,maxy);
end
A = -4*eye((nx-2)*(ny-2),(nx-2)*(ny-2));
b = ones((nx-2)*(ny-2),1).*(-1);
for i=1:(nx-2)*(ny-2)
if mod(i,nx-2) == 1
if i==1
A(1,2) = 1;
A(1,nx-1) = 1;
b(1) = - u0(1,2) - u0(2,1);
else
if i == (ny-3)*(nx-2)+1
A(i,i+1) = 1;
A(i,i-nx+2) = 1;
b(i) = - u0(ny-1,1) - u0(ny,2);
else
A(i,i+1) = 1;
A(i,i-nx+2) = 1;
A(i,i+nx-2) = 1;
b(i) = - u0(floor(i/(nx-2))+2,1);
end
end
else
if mod(i,nx-2) == 0
if i == nx-2
A(i,i-1) = 1;
A(i,i+nx-2) = 1;
b(i) = - u0(1,nx-1) - u0(2,nx);
else
if i == (ny-2)*(nx-2)
A(i,i-1) = 1;
A(i,i-nx+2) = 1;
b(i) = - u0(ny-1,nx) - u0(ny,nx-1);
else
A(i,i-1) = 1;
A(i,i-nx+2) = 1;
A(i,i+nx-2) = 1;
b(i) = - u0(floor(i/(nx-2))+1,nx);
end
end
else
if i>1 && i< nx-2
A(i,i-1) = 1;
A(i,i+nx-2) = 1;
A(i,i+1) = 1;
b(i) = - u0(1,i+1);
else
if i > (ny-3)*(nx-2) && i < (ny-2)*(nx-2) A(i,i-1) = 1;
A(i,i-nx+2) = 1;
A(i,i+1) = 1;
b(i) = - u0(ny,mod(i,(nx-2))+1);
else
A(i,i-1) = 1;
A(i,i+1) = 1;
A(i,i+nx-2) = 1;
A(i,i-nx+2) = 1;
end