I'm trying to use the five-point finite-difference method in Matlab to solve the boundary value problem:

enter image description here

I've spotted a solution:

enter image description here

Using the following code

function Poisson_diko_m(n,niter)

% Solves the Poisson equation
% using the five-point finite-difference method

b=zeros(niter,1);

for iter=1:niter
h=1/n;n2=(n-1)^2;
%fprintf('n =%2.0f\n',n)

% Compute grid
x=zeros(n+1);y=zeros(n+1);
x=h*(0:n);y=h*(0:n);   

% matrix of grid values
u=zeros(n+1,n+1);

% Boundary conditions
u(1,:)=0;
u(n+1,:)=0;
u(:,1)=0;
u(:,n+1)=0;

%Matrix A and identity matrix
A=zeros(n-1,n-1);      
A=sparse((-2/h^2)*diag(ones(n-1,1))+(1/h^2)*diag(ones(n-2,1),1)+(1/h^2)*diag(ones(n-2,1),-1));   
%Full matrix
  AA=sparse(kron(A,eye(n-1))+kron(eye(n-1),A));

%right hand side
%for i=1:n-1
%for j=1:n-1
%ij=(i-1)*(n-1)+j;
%G(ij,1)=fr1(x(i+1),y(j+1));  
%end 
%end
G=reshape(fr1(x(2:n),y(2:n)),(n-1)^2,1);

% Solution of global system      
 C=AA\G;
 % Recover coefficients
  u(2:n,2:n)=(reshape(C,n-1,n-1))';
 %disp('______________________________________________________________________')
 %fprintf('   t           x               u           exact       error')
 %fprintf('\n')
 %disp('______________________________________________________________________')

 E=f1(x,y)';
%for i=1:n+1
%for j=1:n+1
%ij=(i-1)*(n-1)+j;
%E(i,j)=f1(x(i),y(j));  
%end 
%end 

  b(iter,1)=max(max(abs(E-u)));
  if iter==1
      rate=0;
  else
      rate=log(b(iter-1,1)/b(iter,1))/log(2);
  end
     fprintf('n=  %6.0f   Maximum error = %12.3e   Rate =     %12.4f\n',n,b(iter,1),rate);
    disp('______________________________________________________________________')
     n=2*n;
end
%for i=2:niter
%    rate(i)=log(b(i-1)/b(i))/log(2);
%end

%rate
A=(E-u)';

size(A);
x=0:h:1;
y=0:h:1;
[j nn]=size(x);
[l mm]=size(y);
z=reshape(A,mm,nn);
[X,Y]=meshgrid(x,y);
surf(X,Y,z);
zlim([1.e-8 1.e-7]);
ylabel('y'),xlabel('x'),zlabel('Error')
colormap(winter)

f = findobj('Type','surface');
set(f,'FaceLighting','phong');

material shiny
shading interp
light


function f=f1(x,y)
f=x.*(1-y').*cos(pi/2*x).*sin(pi/2*y);


function fr=fr1(x,y)
fr=(-1)*pi*(1-y').*sin(pi/2*x)-pi*x.*cos(pi/2*x).*cos(pi/2*y)-pi^2/2*x.*(1-y').*cos(pi/2*x).*sin(pi/2*y);

I called the function

Poisson_diko_m(2^3,8)

because I wanted approximations for h=1/2^l, l = 3,...,10

But the results seem to be wrong as I expected an upside down cup for the graph and for the error to be much smaller, but instead I got these results:

enter image description here

enter image description here

I'm guessing the mistake must be in how I implemented the two formulae in the far bottom of the Matlab code and not in the rest code. Maybe the multiplications dimensions are wrong but I'm not sure. I tried some dummy x and y to see the results but they seem normal. I do not know where the mistake is. Any help is appreciated.

0 Answers