- Don't use anonymous functions unless absolutely necessary; they are not needed in this problem!
- Vectorize. You don't need to construct A & b in a for loop.
- 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.
- Use exp(), not e^

# Solve the generalized form of the Poisson equation

28 views (last 30 days)

Show older comments

Samuel Martinez
on 20 Jun 2021

Commented: Samuel Martinez
on 20 Jun 2021

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.

##### 0 Comments

### Accepted Answer

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:

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

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!