偏微分方程上机实验报告.doc

  • 格式:doc
  • 大小:72.50 KB
  • 文档页数:6

下载文档原格式

  / 6
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

上机实验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