Solving the Laplace Equation (Dirichlet case) using the Discrete Sine Transform

Contents

You first have to download m-files

Download the m-files sintr.m , isintr.m and put them in the same directory as your other m-files.

Solve the Laplace equation $u_{xx} + u_{yy} = 0$ in a rectangle

We want to find a function $u(x,y)$ for $x\in [0,A]$, $y\in[0,B]$ such that

$$u_{xx}(x,y) + u_{yy}(x,y) = 0 \qquad x\in (0,A)$, $y\in(0,B) \quad\mbox{(PDE)}$$

We prescribe Dirichlet boundary conditions on all four sides of the rectangle:

$$ u(x,0) = 0 , \quad u(x,B) = g(x) \qquad x\in (0,A)$$

$$ u(0,y) = 0 , \quad u(A,y) = 0 \qquad y\in (0,B)$$

At the bottom, left, and right side of the rectangle we want $u=0$. At the top side of the rectangle we want $u(x,B)=g(x)$ with a given function $g(x)$ for $x\in (0,A)$.

The method of Separation of Variables gives the following Solution Formula:

First find the sine series $g(x)=\sum_{k=1}^\infty b_k \sin(k\pi x/A)$. Let $\mu_k:=k\pi/A$ and $C_k:=b_k/\sinh(\mu_k B)$ for $k=1,2,\ldots$. Then the solution $u(x,y)$ is given by

$$u(x,y) = \sum_{k=1}^\infty C_k \sin(\mu_k x)\cdot \sinh(\mu_k y)$$

Example: For a rectangle with lengths $A=3$ and $B=2$ use the given values $g(x)=\frac32-|x-\frac32|$ at the top side of the rectangle.

A=3; B=2;                      % lengths of rectangle sides
fct = @(x) 1.5-abs(x-1.5);     % given function for top side of rectangle y=B
N = 60;                        % use spacing A/N for x
x = A*(1:N-1)'/N;              % only interior nodes            (column vector)
g = fct(x);                    % evaluate g for 1/N,...,(N-1)/N (column vector)

b = sintr(g);                  % take Discrete Sine Transform
kv = (1:N-1)';                 % column vector of k-values
muv = kv*pi/A;                 % column vector of mu-values
C = b./sinh(muv*B);            % coefficients C_k from solution formula

M = 30;
y = B*(0:M)/M;                 % find solution for these y-values (row vector)

W = sinh(muv*y);               % array of sinh(mu_k*y) for k=1...N-1, for all y-values
                               % (column vector)*(row vector) gives array
U = isintr(diag(C)*W);         % multiply 1st row by C(1), 2nd row by C(2), ...
                               % then take Inverse Sine Transform of each column

surf(x,y,U');                  % make surface plot of solution
xlabel('x'); ylabel('y'); title('u(x,y)')
view(-15,30);


Published with MATLAB® R2018a