偏微分方程∂^2*u/∂x^2+∂^2*u/∂y^2=2(x^2+y^2)的数值解法..定解条件 u（0,y）=u（x,0）=0；u（1,y）=y^2;u(x,1)=x^2；我知道精确解是u(x,y)=(xy)^2.但是我想知道该方程的数值解法,和迭代方法.

偏微分方程∂^2*u/∂x^2+∂^2*u/∂y^2=2(x^2+y^2)的数值解法..定解条件 u（0,y）=u（x,0）=0；u（1,y）=y^2;u(x,1)=x^2；我知道精确解是u(x,y)=(xy)^2.但是我想知道该方程的数值解法,和迭代方法.
偏微分方程∂^2*u/∂x^2+∂^2*u/∂y^2=2(x^2+y^2)的数值解法..
定解条件 u（0,y）=u（x,0）=0；u（1,y）=y^2;u(x,1)=x^2；我知道精确解是u(x,y)=(xy)^2.
但是我想知道该方程的数值解法,和迭代方法.比如求 （x,y）=（0.2,0.5）时u的值.不用精确解.

偏微分方程∂^2*u/∂x^2+∂^2*u/∂y^2=2(x^2+y^2)的数值解法..定解条件 u（0,y）=u（x,0）=0；u（1,y）=y^2;u(x,1)=x^2；我知道精确解是u(x,y)=(xy)^2.但是我想知道该方程的数值解法,和迭代方法.

源程序：

function [u,x,y] =Helmholtz(f,g,bx0,bxf,by0,byf,D,Mx,My,MinErr,MaxIter)

%解方程： u_xx + u_yy + g(x,y)u = f(x,y)

% 自变量取值区域 D = [x0,xf,y0,yf] = {(x,y) |x0 <= x <= xf, y0<= y <= yf}

% 边界条件

% u(x0,y) = bx0(y), u(xf,y) = bxf(y)

% u(x,y0) = by0(x), u(x,yf) = byf(x)

% x轴均分为Mx段

% y轴均分为My段

% tol 误差因子

% MaxIter: 最大迭代次数

%[u,x,y]:方程u(x,y)在(x,y)点的函数值

x0 = D(1); xf = D(2); y0 = D(3); yf = D(4);

dx = (xf - x0)/Mx; x = x0 + [0:Mx]*dx;%构造内点数组

dy = (yf - y0)/My; y = y0 + [0:My]'*dy;

Mx1 = Mx + 1; My1 = My + 1;

%边界条件

for m = 1:My1

    u([1Mx1],m)=[bx0(y(m)) bxf(y(m))]; %左右边界

end

for n = 1:Mx1

    u(n,[1 My1]) =[by0(x(n)); byf(x(n))];%上下边界

end

%边界平均值作迭代初值

sum_of_bv = sum(sum([u([1 Mx1],2:My) u(2:Mx,[1 My1])']));

u(2:Mx,2:My) = sum_of_bv/(2*(Mx + My - 2));

for i = 1:Mx

    for j = 1:My

         F(i,j) =f(x(i),y(j)); G(i,j) = g(x(i),y(j));

    end

end

dx2 = dx*dx; dy2 = dy*dy; dxy2 = 2*(dx2 + dy2);

rx = dx2/dxy2; ry = dy2/dxy2; rxy = rx*dy2;

for itr = 1:MaxIter

    for i = 2:Mx

        for j =2:My

            u(i,j)= ry*(u(i+1,j)+u(i-1,j)) + rx*(u(i,j+1)+u(i,j-1))+ rxy*(G(i,j)*u(i,j)- F(i,j));%迭代公式

        end

    end

    if itr > 1& max(max(abs(u - u0))) < MinErr%循环结束条件

        break;

    end

    u0 = u;

end

u=u';

在MATLAB中编写脚本文件：

f = inline('2*x^2+2*y^2','x','y');

g = inline('0','x','y');

x0 = 0; xf = 1; y0 = 0; yf = 1;%自变量取值范围

Mx = 50;My = 30;%等分段数

bx0 = inline('0','y'); %边界条件

bxf = inline('y^2','y');

by0 = inline('0','x');

byf = inline('x^2','x');

D = [x0 xf y0 yf]; MaxIter = 100; MinErr = 1e-4;

[U,x,y] =Helmholtz(f,g,bx0,bxf,by0,byf,D,Mx,My,MinErr,MaxIter);

clf, mesh(U)

xlabel('x')

ylabel('y')

zlabel('U')