Solve the generalized form of the Poisson equation

Hello, I am trying to solve the following problem in a rectangle with Dirichlet conditions at the boundary. I have the following implementation for this problem:

n =25; 
dx = 1/(n-1); 
x= 0:dx:1;
y= x;
[X,Y] = ndgrid(x,y);
 

isDirichlet = (X==0) | (X==1) | (Y==0) | (Y==1);
 
u = @(x,y) 10*x.*y.*(1-x).*(1-y).*exp(x.^(4.5)); 
k = @(x,y) cos(x);
f = @(x,y) e^(x.^4.5).*(-45 *(-x.^5.5 + x.^4.5 - 0.444444*x + 0.222222).* (y - 1).* y.* sin(x)...
     + 247.5* (-1.36364*x.^4.5 + x.^3.5 - 0.818182*x.^9 + 0.818182*x.^8 - 0.0808081).*(y - 1).*y.*cos(x)...
     - 20*(x - 1).*x.*cos(x)); 
g = @(x,y) 0; 

A = spalloc(n^2,n^2,5);
b = zeros(n^2,1);
 for i=1:n^2
l(i)=isDirichlet(i);
end
 

for i = 1:n^2
   if isDirichlet(i)
       A(i,i) = 1;
       b(i) = g(X(i),Y(i));
       continue
   end
   
   [row,col] = ind2sub([n,n],i); 
   L = sub2ind([n,n],row-1,col); 
   R = sub2ind([n,n],row+1,col); 
   U = sub2ind([n,n],row,col+1); 
   D = sub2ind([n,n],row,col-1);
   
   kl = k((X(i)+X(L))/2,Y(i));
   kr = k((X(i)+X(R))/2,Y(i));
   ku = k(X(i),(Y(i)+Y(U))/2);
   kd = k(X(i),(Y(i)+Y(D))/2);

   A(i,i) = kl+kr+ku+kd;
   A(i,L) = -kl;
   A(i,R) = -kr;
   A(i,U) = -ku;
   A(i,D) = -kd;
   
   b(i) = f(X(i),Y(i))*dx*dx;
end

uh=A\b;
Uh = reshape(uh,[n,n]);
surf(Uh); hold on
surf(u(X,Y)+1);
max(max(abs(Uh-u(X,Y))))

With the above it gives me the following error. With n = 50 I have the following error::

error: sub2ind: index (51,_): out of bound 50 
error: called from 
prueba at line 33 column 6

But with n = 51 it still works, I don't know what this error could be. Also, I have to solve this for 500, 1000, 5000 and 10000 points, but for that I would need a lot of memory, so I guess you can take advantage of the matrix structure in some way to achieve these results. I appreciate if you could help me to fix that error and see how I can implement for 500, 1000, 5000 and 10000. First of all, thank you.

Answer 1

[Joel Lynch]

I don't get an error when running your code with MATLAB version 2019b, but the bigger problem you have is that your method will not scale well. A few tips off the top:

1. Don't use anonymous functions unless absolutely necessary; they are not needed in this problem!
2. Vectorize. You don't need to construct A & b in a for loop.
3. Use spdiags. For a 5-point stencil there are 5 non-zero elements in each row (forming 5 diagonals), and you can use spdiags to generate your sparse matrix from an Nx5 matrix, rather than using spalloc and filling entries manually.
4. Use exp(), not e^

I put together the following solution, which runs values 50, 100, and 500 in 0.027, 0.11, ands 1.56 seconds respectively on my computer. Your code runs the same values in 0.074, 0.31, and 405 seconds. That's >250 times faster for N=500!

close all
clear all
clc

% Grid Definition
nx = 50;
N = nx*nx; 
x = linspace(0,1,nx);
y = x;
[X,Y] = ndgrid(x,y);
dx = x(2)-x(1);

% Compute RHS with no BC, assuming f(X,Y)
b = -1.36364*X.^4.5 +  X.^3.5 - 0.818182*X.^9 + 0.818182*X.^8 - 0.0808081;
b = b * 247.5 .* (Y - 1).*Y.*cos(X);
b = b -45 *(-X.^5.5 + X.^4.5 - 0.444444*X + 0.222222).*(Y - 1).* Y.* sin(X);
b = b -20 *(X - 1).*X.*cos(X) ; 
b = exp( X.^4.5 ) .* b;
b = b.*dx.*dx;

% Overwrite BC's on the RHS using vectorized assignment
isDirichlet = (X==0) | (X==1) | (Y==0) | (Y==1);
b(isDirichlet) = 0.0;

% Compute A Matrix diagonals Ad(1:N,1:5) assuming normal stencils
Ad = zeros(N,5);

% Compute off-center points, note that there is no y dependence in
% k(x,y)=cos(x)
Ad(:,1) = -cos(X(:)); % Down = -kd = -cos(X)
Ad(:,2) = -cos(X(:)-dx*0.5); % Left = -kl = -cos(X-dx/2)
Ad(:,4) = -cos(X(:)+dx*0.5); % Right= -kl = -cos(X+dx/2)
Ad(:,5) = Ad(:,1); % Up = -ku = -cos(X) = Down

% Use these to get center by summing the negatives
Ad(:,3) = -sum( Ad(:,[1,2,4,5]), 2 ); % Center

% Overwrite BC's on the A Matrix diagonals
idx = find(isDirichlet(:)); % array of inidices of BC points
Ad(idx,:) = 0.0; % Zero all elements of rows defined by idx
Ad(idx,3) = 1.0; % Make center node=1

% Create A Matrix from Ad, assigning columns to diagonals
A = spdiags(Ad,[-nx,-1,0,1,nx],N,N)'; % ADDED TRANSPOSE

% Solve & Reshape
uh=A\b(:);
uh = reshape(uh,[nx,nx]);

% Compute Exact solution and error
u_exact = 10*X.*Y.*(1-X).*(1-Y).*exp(X.^(4.5)); 
u_error = abs( (uh - u_exact)./u_exact );

% Reshape and plot solution on two subplots
subplot(1,3,1); view(0,90); hold on; colorbar;
surf(x,y,uh,'edgecolor','none');
xlabel('x'); ylabel('y')
title('Computed')
shading interp
axis square

% Add exact solution
subplot(1,3,2); view(0,90); hold on; colorbar; 
surf(x,y,u_exact,'edgecolor','none');
xlabel('x'); 
title('Exact')
shading interp
axis square

% add Error
subplot(1,3,3); view(0,90); hold on; colorbar;
surf(x,y,u_error,'edgecolor','none');
xlabel('x'); 
title('Error')
shading interp
axis square

[Samuel Martinez]

You are right, in matlab it does not fail for n = 50, which is strange, I was working in octave and that was where it had the failure. I put the implementation the way I show in my code because I am going to change f and k for example to k = 1 + x ^ 2 + y ^ 2 and choose f in such a way that the exact solution is the same as for the previous implementation. I notice that if your code works much faster, but is there a way to adapt it for the general problem? Thank you very much for your answers in advance.

[Joel Lynch]

You should be able to drop in your original anonymous function calls when defining b and k. As long as you are still vectorizing and not calling those functions repeatedly in a loop, there shouldn't be much of a performance hit for using them.

[Joel Lynch]

As an example: you could alter the Ad section to be:

k = @(x,y) 1 + x .^ 2 + y .^ 2; 
Ad(:,1) = -k( X(:), Y(:)-dx/2 );   % Down 
Ad(:,2) = -k( X(:)-dx/2, Y(:) );   % left
Ad(:,4) = -k( X(:)+dx/2, Y(:) );   % right 
Ad(:,5) = -k( X(:), Y(:)+dx/2 );   % Up 

Note that you should probably give k a name like "k_func" as single letters are usually indicies.

[Samuel Martinez]

I think there is an error, it does not put uh = 0 on the border, for the simple case of n = 3 this phenomenon happens

[Joel Lynch]

Thanks for catching that; the spdiags command works a little different than I expected. Rather than construct each row of A from each row of Ad, it constructs each column of A from each row of Ad. The solution is simply to transpose the result of the spdiags function. I added that to the solution.

[Samuel Martinez]

I think there is something wrong, I can't find it, in theory the error should decrease by 1/4 if I double the mesh, which does not happen here. Or I'm wrong?

The question that I had has already been solved and thank you very much, these questions that I am asking you are extra and I thank you for answering.

[Joel Lynch]

To clarify, the line A = spdiags(Ad,[-nx,-1,0,1,nx],N,N)'; has an added apostrophe on the end, which transposes A.

I can't say if the exact solution is correct and implemented properly, but that could be a source of confusion if it is not correct. If indeed the computed solution is not converging like a second-order scheme, I would check that your stencil is being applied correctly, as well as the "f" function, which is long and could contain an error.

Good luck!

[Samuel Martinez]

Thank you. Greetings!