Search code examples
matlabmathpde

Finite difference method for elliptic equations

If R is the region of the plane (0,1) x (0,2), Let L be the 2-d Laplace operator and consider the Poisson equation Lu = 4 on R. One solution is the function v(x, y) = (x-y)^2. Let g be the restriction of v to the boundary of R. To obtain it, we set h = k = 1/2, and m=3, n =5.

Here is my code: but it doesn't work.

function w=poisson(xl,xr,yb,yt,M,N)
f=@(x,y) 0; % define input function data
g1=@(x) (x^2-2*x+1); % define boundary values
g2=@(x) (x^2-4*x+4); % Example 8.8 is shown
g3=@(y) y^2;
g4=@(y) y^2;
m=M+1;n=N+1; mn=m*n;
h=(xr-xl)/M;h2=h^2;k=(yt-yb)/N;k2=k^2;
x=xl+(0:M)*h; % set mesh values
y=yb+(0:N)*k;
A=zeros(mn,mn);b=zeros(mn,1);
for i=2:m-1 % interior points
  for j=2:n-1
    A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;A(i+(j-1)*m,i+1+(j-1)*m)=1/h2;
    A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2;
    A(i+(j-1)*m,i+(j-2)*m)=1/k2;A(i+(j-1)*m,i+j*m)=1/k2;
    b(i+(j-1)*m)=f(x(i),y(j));
  end
end
for i=1:m % bottom and top boundary points
  j=1;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g1(x(i));
  j=n;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g2(x(i));
end
for j=2:n-1 % left and right boundary points
  i=1;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g3(y(j));
  i=m;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g4(y(j));
end
v=A\b; % solve for solution in v labeling
w=reshape(v(1:mn),m,n);  %translate from v to w
mesh(x,y,w')

enter image description here

Solution

  • As per your problem description, I changed the following and got an answer:

    h=1/2;