Solve the generalized form of the Poisson equation

13 views (last 30 days)
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.

Accepted Answer

Joel Lynch
Joel Lynch on 20 Jun 2021
Edited: Joel Lynch on 20 Jun 2021
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
  8 Comments
Joel Lynch
Joel Lynch on 20 Jun 2021
Edited: Joel Lynch on 20 Jun 2021
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!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!