# Solve the generalized form of the Poisson equation

21 visualizaciones (últimos 30 días)
Samuel Martinez el 20 de Jun. de 2021
Comentada: Samuel Martinez el 20 de Jun. de 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 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

Joel Lynch el 20 de Jun. de 2021
Editada: Joel Lynch el 20 de Jun. de 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
% 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)
% Use these to get center by summing the negatives
% 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
% 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')
axis square
subplot(1,3,2); view(0,90); hold on; colorbar;
surf(x,y,u_exact,'edgecolor','none');
xlabel('x');
title('Exact')
axis square
subplot(1,3,3); view(0,90); hold on; colorbar;
surf(x,y,u_error,'edgecolor','none');
xlabel('x');
title('Error')
axis square
##### 8 comentariosMostrar 6 comentarios más antiguosOcultar 6 comentarios más antiguos
Joel Lynch el 20 de Jun. de 2021
Editada: Joel Lynch el 20 de Jun. de 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!
Samuel Martinez el 20 de Jun. de 2021
Thank you. Greetings!

Iniciar sesión para comentar.

### Categorías

Más información sobre Geometry and Mesh en Help Center y File Exchange.

### Community Treasure Hunt

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

Start Hunting!

Translated by