Question

我有个偏微分方程的求解问题。这题目应该是用二阶差分方式解得吧。。但是就不知道怎么编程。麻烦您了

在单位正方形上有偏微分方程D^2U/Dx^2+D^2U/Dy^2=2(x^2+y^2)，定解条件是u（0，y）=u(x,0)=0;u(1,y)=y^2,u(x,1)=x^2.
分别取N=5,10,15把单位正方形划分成N*N个小格子，并把格子线的交点编号（i,j），
记u(x(i),x(j))=u(i,j).用熟知微分公式建立关于u（i,j）的线性代数方程组，并求解。
该方程的精确解是u（x，y）=（xy）^2，将精确解和近似解画在同一张三维图片上。
题目就是这样了，麻烦您帮帮忙啊。。我实在解不出来。。拜托了。。万分万分感激啊..我的QQ号是61972352。。邮箱是aprilnan@126.com。。

Answer 1

function varargout=liu(varargin)
N=5;%改这个N=10 15
a=0;b=1;c=0;d=1;h=1/(N-1);
f=inline('-2*(x^2+y^2)','x','y');
g1x=inline('0');
g2x=inline('x^2');
g1y=inline('0');
g2y=inline('y^2');
chfenmethed(f,g1x,g2x,g1y,g2y,a,b,c,d,h);
function [X,Y,Z]=chfenmethed(f,g1x,g2x,g1y,g2y,a,b,c,d,h)
%求解下问题
%u_xx+u_yy=-f(x,y) x,y 在区域内 x in a<x<b,c<y<d
%u=g(x,y) x,y在边界上
%u=g(a,y)=g1y u=g(b,y)=g2y c=<y<=d
%u=g(x,c)=g1x u=g(x,d)=g2x
a<x<b
%h离散x方向的步长
%h离散y方向的步长
N=10000;
x=a:h:b;
y=c:h:d;
m=length(x);
n=length(y);
ee=0.00001;
[Y,X]=meshgrid(y,x);
Uliu=f0(X,Y);
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

s1=U(i+1,j)+U(i,j+1)+U(i-1,j)+U(i,j-1);

s2=U(i+1,j+1)+U(i-1,j+1)+U(i+1,j-1)+U(i-1,j-1);

s3=feval(f,x(i+1),y(j))+feval(f,x(i-1),y(j))...

+feval(f,x(i),y(j+1))+feval(f,x(i),y(j-1));

U(i,j)=(4*s1+s2-h^2/2*s3-4*h^2*feval(f,x(i),y(j)))/20;

end

end
if max(max(abs(U0-U)))<ee

break
end
end
%-----------------------------------------------------
mesh(X,Y,U)
hold on
plot3(X,Y,Uliu,'ro')
function z=f0(x,y)
%精确解函数
z=(x.*y).^2;