用五点有限差分格式求解椭圆型方程(偏微分方程) 程序2

用五点有限差分格式求解椭圆型方程（偏微分方程）程序2

2010-04-29 10:33

function varargout=liu(varargin)

a=0;b=2;c=0;d=1;h1=1/16;h2=1/16;

f=inline('(pi^2-1)*exp(x)*sin(pi*y)','x','y');

g1x=inline('0');

g2x=inline('0');

g1y=inline('sin(pi*y)');

g2y=inline('exp(2)*sin(pi*y)');

[X,Y,Z]=chfenmethed(f,g1x,g2x,g1y,g2y,a,b,c,d,h1,h2);

mesh(X,Y,Z);

shading flat;

xlabel('X','FontSize',14);

ylabel('Y','FontSize',14);

zlabel('error','FontSize',14);

title('误差图');

function [X,T,Z]=chfenmethed(f,g1x,g2x,g1y,g2y,a,b,c,d,h1,h2) %求解下问题

%-(u_xx+u_yy)=f(x,y) x,y 在区域内x in a

%u=g(x,y) x,y在边界上

%u=g(a,y)=g1y u=g(b,y)=g2y c=

%u=g(x,c)=g1x u=g(x,d)=g2x a

%h1离散x方向的步长

%h2离散y方向的步长

N=10000;

x=a:h1:b;

y=c:h2:d;

m=length(x);

n=length(y);

ee=0.00001;

[X,T]=meshgrid(x,y);

Z=zeros(n,m);

U=zeros(n,m);

for i=2:m-1

U(1,i)=feval(g1x,x(i));

U(n,i)=feval(g2x,x(i));

end

for j=1:n

U(j,1)=feval(g1y,y(j));

U(j,m)=feval(g2y,y(j));

end

%while true

%下为高斯赛德尔迭代法

%---------------------------------------------------------------------- for k=1:N

U0=U;

for i=2:n-1

for j=2:m-1

U(i,j)=(feval(f,x(j),y(i))+1/h2^2*U(i-1,j)+1/h1^2*U(i,j-1)...

+1/h1^2*U(i,j+1)+1/h2^2*U(i+1,j))/(2/h1^2+2/h2^2);

end

end

if max(max(abs(U0-U)))

break

end

end

%-----------------------------------------------------

for i=2:n-1

for j=2:m-1

qu=f0(x(j),y(i));

Z(i,j)=abs(qu-U(i,j));

end

end

function z=f0(x,y)

%精确解函数

z=exp(x)*sin(pi*y);