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

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

2010-04-20 22:31

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,Y,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.0001;

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

Z=zeros(m,n);

U=zeros(m,n);

for i=2:m-1

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

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

end

for j=1:n

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

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

end

%while true

%下为雅克比迭代法

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

U0=U;

for i=2:m-1

for j=2:n-1

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

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

end

end

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

break

end

end

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

for i=2:m-1

for j=2:n-1

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

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

end

end

function z=f0(x,y)

%精确解函数

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