% Elliptic Partial Differential Equations - Iterative 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
% Solving this for  T_{i,j}we obtain:
%  T_{i,j} =\frac{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)}
% To solve this system using iteration, we set initial (guess) values for all of the temperatures and then go through the mesh points one-by-one to update to new values using the equation above. We then repeat the process until we've reached convergence.

clear all
Q=-1e5; %volumetric heating
cond=100; %thermal conductivity
width=1; height=1;
%set boundary values
w_top=0; w_left=0; w_right=0; w_bottom=0;
%set number of divisions in both directions
ndiv=100; n=ndiv+1; mdiv=100; m=mdiv+1;
centeri=round(ndiv/2+1); centerj=round(mdiv/2+1);
%calculate width of divisions in both directions
h=width/mdiv; s=height/ndiv;
x=linspace(0,width,n); y=linspace(0,height,m);
w=zeros(n,m,'double');
for i=1:n
    w(i,1)=w_left; w(i,m)=w_right;
end
for j=1:m
    w(1,j)=w_top; w(n,j)=w_bottom;
end
%set initial values for iteration
for i=2:n-1
    for j=2:m-1
        w(i,j)=0;
    end
end
wpeakold=0;
numiters=10000;
relativechange=zeros(numiters,1);
for k=1:numiters
    for i=2:n-1
        for j=2:m-1
            w(i,j)=(w(i,j+1) + w(i,j-1) + h^2/s^2*(w(i-1,j) + w(i+1,j))-Q/cond*h^2)/(1+h^2/s^2)/2;
        end
    end
    %track relative change in peak value of w to track convergence
    newwpeak=w(centeri,centerj);
    relativechange(k)=(abs(wpeakold-newwpeak)/newwpeak);
    wpeakold=newwpeak;
end
Maximum_value_iteration=newwpeak
semilogy(1:numiters,relativechange, 'LineWidth',2)
xlabel('iteration number'); ylabel('relative change in peak value');
set(gca,'FontSize',10)
mesh(x,y,w')
xlabel('X'); ylabel('Y');
contour(x,y,w','ShowText','on');
xlabel('X');ylabel('Y')

% Note that I needed about 10,000 iterations to get reasonably close to the previously determined solution for the peak temperature. If I make changes to the model, I might require more (or fewer) iterations.
% We can speed up the convergence by using successive overrelaxation. The idea here is that as we iterate to the solution we artificially increase the increment from one iteration to the next. That is, if w_{i,j} is the value of the temperature prior to the current iteration and {wnew}_{i,j} is the updated value of the temperature using the iteration scheme derived previously, then we use the following as our new temperature:
% T_{i,j} =w_{i,j} +L*\left({wnew}_{i,j} -w_{i,j} \right)
% In the code below we use L=1.5, implying that as we iterate we move about 50% farther than the original algorithm would suggest. The code is as follows:

clear all
Q=-1e5; %volumetric heating
cond=100; %thermal conductivity
width=1; height=1;
%set boundary values
w_top=0; w_left=0; w_right=0; w_bottom=0;
%set number of divisions in both directions
ndiv=100; n=ndiv+1; mdiv=100; m=mdiv+1;
centeri=round(ndiv/2+1); centerj=round(mdiv/2+1);
%calculate width of divisions in both directions
h=width/mdiv; s=height/ndiv;
x=linspace(0,width,n); y=linspace(0,height,m);
L=1.5;      %overrelaxation value
w=zeros(n,m,'double');
w1=zeros(n,m,'double');
for i=1:n
    w(i,1)=w_left; w(i,m)=w_right;
end
for j=1:m
    w(1,j)=w_top; w(n,j)=w_bottom;
end
%set initial values for iteration
for i=2:n-1
    for j=2:m-1
        w(i,j)=0;
    end
end
wpeakold=0;
numiters=4000;
relativechange=zeros(numiters,1);
for k=1:numiters
    for i=2:n-1
        for j=2:m-1
            w1(i,j)=(w(i,j+1) + w(i,j-1) + h^2/s^2*(w(i-1,j) + w(i+1,j))-Q/cond*h^2)/(1+h^2/s^2)/2;
            w(i,j)=L*w1(i,j)+(1-L)*w(i,j);
        end
    end
    %track relative change in peak value of w to track convergence
    newwpeak=w(centeri,centerj);
    relativechange(k)=(abs(wpeakold-newwpeak)/newwpeak);
    wpeakold=newwpeak;
end
Maximum_value_overrelaxation=newwpeak
semilogy(1:numiters,relativechange, 'LineWidth',2)
xlabel('iteration number'); ylabel('relative change in peak value');
set(gca,'FontSize',10)
mesh(x,y,w')
xlabel('X'); ylabel('Y');
contour(x,y,w','ShowText','on');
xlabel('X');ylabel('Y')

% Notice that I used 4,000 iterations in the code using overrelaxation, as opposed to 10,000, and got slightly better results.
% © James Blanchard, University of Wisconsin-Madison, blanchard@engr.wisc.edu