% Elliptic Partial Differential Equations - Direct Approach
% Elliptic partial differential equations can be thought of as two-dimensional boundary value problems. The classic problem is heat conduction in a square plate, where boundary conditions are provided on the 4 sides of a plate and we seek the temperature distribution in the interior of the plate. Other examples of elliptic PDE’s include diffusion, electric fields, and torsion.  We will deal, initially, with steady-state heat conduction problems in two dimensions.  In these problems, we will seek solutions for the temperature field in some 2-D region with specified boundary conditions.  The governing equation for these problems (in Cartesian coordinates) is:
%  \frac{\partial^2 T}{\partial x^2 }+\frac{\partial^2 T}{\partial y^2 }+\frac{Q}{k}=0
% where T is the temperature, Q is the volumetric heating, and k is the thermal conductivity of the solid.
% As a model problem, we’ll consider a rectangular plate of width a and height b, so our region of interest is the rectangle defined by:  0<x<a and 0<y<b.  For boundary conditions, we will assume we know the temperature along the entire boundary of the plate is zero.  These are known as Dirichlet boundary conditions. This equation can be solved using separation of variables, giving the solution for the maximum temperature (at the center of the plate) as:
%  \begin{array}{l}
% T_{\max } =\left(\frac{1}{8}+\frac{1}{2}\sum_{k=1}^{\infty } \frac{{\left(-1\right)}^k }{\lambda_k^3 \;\cosh \left(\lambda_k \right)}\right)\frac{Q}{k}\\
% \lambda_k =\frac{1}{2}\left(2k-1\right)\pi
% \end{array}
% For comparison to the numerical results below, this evaluates to:

clear all
Q=1e5; cond=100;
numterms=10;
kvec=1:numterms;
lam_k=1/2*(2*kvec-1)*pi;
Tmax_analytical=Q/cond*(1/8+1/2*sum((-1).^kvec./(lam_k.^3.*cosh(lam_k))))

% To solve this problem numerically, we can use a technique called the finite difference method. We divide the rectangle into a grid and assign a temperature to each grid point.  This is the 2-D equivalent of dividing a 1-D region into intervals to create mesh points.  We will label each temperature using the notation T_{i,j}, where i denotes the row number and j denotes the column number.  We will assume that the grid spacing in the x-direction is given by h and the spacing in the y-direction is s.  Given this model, we can approximate the second partial derivative of the temperature with respect to x as:
%  \frac{\partial^2 T}{\partial x^2 }\approx \frac{T_{i,j+1} -2T_{i,j} +T_{i,j-1} }{h^2 }
% and the second partial derivative of  the temperature with respect to y as:
%  \frac{\partial^2 T}{\partial y^2 }\approx \frac{T_{i+1,j} -2T_{i,j} +T_{i-1,j} }{s^2 }
% If we substitute these into our PDE and set Q=0, then we obtain:
%   \frac{T_{i,j+1} -2T_{i,j} +T_{i,j-1} }{h^2 }+\frac{T_{i+1,j} -2T_{i,j} +T_{i-1,j} }{s^2 }+\frac{Q}{k}=0
% To solve this problem using the direct approach, we first rewrite this equation such that the unknown quantities are on the left and the known quantities are on the right:
% T_{i,j+1} +T_{i,j-1} +\frac{h^2 }{s^2 }\left(T_{i+1,j} +T_{i-1,j} \right)-2\left(1+\frac{h^2 }{s^2 }\right)T_{i,j} =\frac{-Qh^2 }{k}
% We can then write this equation at all the internal mesh points, account for the boundary conditions, and solve the resulting set of algebraic equations. This is achieved in the code below.

clear all
width=1; height=1;
%set up number of divisions and mesh points in each direction
ndiv=20; mdiv=20;
Q=1e5; cond=100;
m=mdiv+1;n=ndiv+1;
h=width/mdiv;s=height/ndiv;
x=(0:mdiv)*h; y=(0:ndiv)*s;
A=zeros(m*n,m*n);b=zeros(m*n,1);
for i=2:m-1 % interior points
    for j=2:n-1
        A(i+(j-1)*m,i-1+(j-1)*m)=1/h^2;
        A(i+(j-1)*m,i+1+(j-1)*m)=1/h^2;
        A(i+(j-1)*m,i+(j-1)*m)=-2/h^2-2/s^2;
        A(i+(j-1)*m,i+(j-2)*m)=1/s^2;
        A(i+(j-1)*m,i+j*m)=1/s^2;
        b(i+(j-1)*m)=-Q/cond;
    end
end
for i=1:m % bottom and top boundary points as well as diagonal of A
    A(i,i)=1;
    b(i)=0;
    A(i+(n-1)*m,i+(n-1)*m)=1;
    b(i+(n-1)*m)=0;
end
for j=2:n-1 % left and right boundary points
    A(1+(j-1)*m,1+(j-1)*m)=1;
    b(1+(j-1)*m)=0;
    A(m+(j-1)*m,m+(j-1)*m)=1;
    b(m+(j-1)*m)=0;
end
v=A\b; % solve for solution in v labeling
w=reshape(v(1:m*n),m,n);  %translate from v to w
mesh(x,y,w')
xlabel('x');ylabel('y');
Tmax_direct=max(max(w))

%
% Note that the discrepancy between the analytical and numerical solutions will decrease as we divide up our rectangle more finely, but the solution will take longer and use more memory. We can improve on both of these fronts by using sparse matrices. This is particularly helpful for these problems because the vast majority of the elements in A are zeros. Specifically, for a mesh that is 20 divisions by 20 divisions, A has 194,481 ({21}^4) total elements and 1,885 (less than 1%) are non-zero.
% The sparse version of our solution is as follows:

clear all
width=1; height=1;
%set up number of divisions and mesh points in each direction
ndiv=200; mdiv=200;
Q=1e5; cond=100;
m=mdiv+1;n=ndiv+1;
h=width/mdiv;s=height/ndiv;
x=(0:mdiv)*h; % set mesh values
y=(0:ndiv)*s;
A=sparse(m*n,m*n);b=zeros(m*n,1);
for i=2:m-1 % interior points
    for j=2:n-1
        A(i+(j-1)*m,i-1+(j-1)*m)=1/h^2;
        A(i+(j-1)*m,i+1+(j-1)*m)=1/h^2;
        A(i+(j-1)*m,i+(j-1)*m)=-2/h^2-2/s^2;
        A(i+(j-1)*m,i+(j-2)*m)=1/s^2;
        A(i+(j-1)*m,i+j*m)=1/s^2;
        b(i+(j-1)*m)=-Q/cond;
    end
end
for i=1:m % bottom and top boundary points as well as diagonal of A
    A(i,i)=1;
    b(i)=0;
    A(i+(n-1)*m,i+(n-1)*m)=1;
    b(i+(n-1)*m)=0;
end
for j=2:n-1 % left and right boundary points
    A(1+(j-1)*m,1+(j-1)*m)=1;
    b(1+(j-1)*m)=0;
    A(m+(j-1)*m,m+(j-1)*m)=1;
    b(m+(j-1)*m)=0;
end
v=A\b; % solve for solution in v labeling
w=reshape(v(1:m*n),m,n);  %translate from v to w
mesh(x,y,w')
xlabel('x');ylabel('y');
Tmax_direct=max(max(w))

% Note that I would not have been able to run a case with 200x200 mesh without sparse, due to memory limitations. The script could be made even more efficient by building the entire A array with sparse.