% Parabolic Partial Differential Equations - Finite Difference Methods
% Parabolic partial differential equations (PDEs) are a class of PDE that is often used to model one-dimensional, time dependent problems. One example is one-dimensional, transient heat conduction. In this file I’ll describe how to solve simple parabolic PDEs using finite difference methods.
% As a model problem, I’ll solve a 1-D transient conduction problem. The problem models a surface heat flux (q”) on an infinite slab with thermal conductivity k, density  \rho, and heat capacity c_p. The slab is of thickness L and the back surface is held at a temperature of 0 degrees for all time. The entire slab is initially at a temperature of 0 degrees and then the surface heat is turned on and held constant. The slab temperatures will rise until a steady state is reached. The equation we’re solving is
% \begin{array}{l}
% \rho c_p \frac{\partial T}{\partial t}=k\frac{\partial^2 T}{\partial x^2 }\\
% T\left(x,t=0\right)=0\\
% -k\frac{\partial T}{\partial x}=\textrm{q"}\;\textrm{at}\;x=0\\
% T=0\;\textrm{at}\;x=L
% \end{array}
% For comparison to our future numerical calculations, the solution for the surface temperature in cases where L is large is as follows:
% T_{\textrm{surface}} =\frac{2\textrm{q"}}{k}\sqrt{\frac{\textrm{kt}}{\rho c_p \pi }}
% To solve with finite difference methods, we first approximate the derivatives with discrete approximations. We divide our object, spatially, into small divisions and approximate the derivatives at the mesh points (the points separating the divisions). We also divide the time into discrete time steps and only solve the system at these discrete points in time. Hence, we will refer to temperatures as T_i^j where i represents the spatial mesh point and j represents the time step. This allows us to represent an approximate time derivative as
% \frac{\partial T}{\partial t}\approx \frac{T_i^{j+1} -T_i^j }{\Delta t}
% where \Delta t is the time step. Similarly, we can write the spatial term as
% \frac{\partial^2 T}{\partial x^2 }\approx \frac{T_{i+1}^j -2T_i^j +T_{i-1}^j }{h^2 }
% where h is the width of our spatial divisions. Note that we've written the spatial derivative at the beginning of the time step. As we'll see, this leads to some issues and we will modify it later. Note that this approach produces an explicit formulation of the problem. Using these approximations, we can rewrite our original differential equation as
% \left(\frac{\rho c_p }{k}\right)\frac{T_i^{j+1} -T_i^j }{\Delta t}\approx \frac{T_{i+1}^j -2T_i^j +T_{i-1}^j }{h^2 }
% Given that only one term here is written at time step j+1, we can solve for that as follows:
% \begin{array}{l}
% T_i^{j+1} \approx T_i^j \left(1-2\alpha \right)+\alpha \left(T_{i+1}^j +T_{i-1}^j \right)\\
% \alpha =\frac{k\Delta t}{\rho c_p h^2 }
% \end{array}
% We can now start with our initial temperatures at every spatial step and use the above algorithm to calculate a new temperate for every mesh point after the first time step. Then we can use those to calculate the spatial temperature distribution after the following time step. This process can continue for as long as we'd like.
% The drawback to the explicit formulation is that this solution is unstable for certain values of \alpha and the error will grow very rapidly. Specifically, we can only obtain stable solutions if \alpha \le \frac{1}{2}.
% To address the boundary conditions, the right side is straightforward because the temperature of that end starts at T=0 and does not change. However, we must address the derivative in the boundary condition at x=0. It is customary to introduce the concept of a "ghost point," which is an imaginary mesh point a distance h outside of the surface. We will call this ghost point T_0^j. In this case we can approximate the spatial derivative at the surface as:
% \frac{\partial T}{\partial x}\left(x=0\right)\approx \frac{T_2^j -T_0^j }{2h}
% We can also rewrite our PDE approximation at the surface with
% T_1^{j+1} \approx T_1^j \left(1-2\alpha \right)+\alpha \left(T_2^j +T_0^j \right)
% We can then combine these two equations to eliminate the ghost point. So our boundary condition is written as
% \begin{array}{l}
% -k\frac{T_2^j -T_0^j }{2h}=\mathrm{q"}\\
% T_0^j =T_2^j +\frac{2h\mathrm{q"}}{k}
% \end{array}
% T_1^{j+1} \approx T_1^j \left(1-2\alpha \right)+2\alpha \left(T_2^j +\frac{h\mathrm{q"}}{k}\right)
% The code would be as follows:

clear all
q=1e5;
k=100; rho=1000; cp=500;
diffus=k/rho/cp;
length=1; 
ndiv=100; %number of spatial divisions
h=length/ndiv; %mesh spacing
endtime=300;
ntimes=1500;
dt=endtime/ntimes;
alpha=diffus*dt/h^2
tboundary=0; %boundary temperature (C)
time=zeros(ntimes+1);
temperature=zeros(ndiv+1);
tsurf=zeros(ntimes+1);
told=zeros(ndiv+1); %fill temperature vector with zeros
for i=1:ntimes
    time(i+1)=i*dt;
    temperature(1)=told(1)*(1-2*alpha)+2*alpha*(told(2)+h*q/k);
    temperature(ndiv+1)=tboundary;
    for j=2:ndiv
        temperature(j)=(1-2*alpha)*told(j)+alpha*(told(j-1)+told(j+1));
    end
    told=temperature;
    tsurf(i+1)=temperature(1);
end
tanalytic=2*q/k*sqrt(time*k/rho/cp/pi);
plot(time,tsurf,time,tanalytic,'o')

% If you want to see what happens with an explicit formulation if \alpha >\frac{1}{2}, try using fewer time steps.
% To avoid these stability issues, we can swith to an implicit formulation by writing our spatial derivative at the end of the time step.
% To implement this explicit formulation, we write the following:
% T_i^{j+1} -T_i^j \approx \alpha \left(T_{i+1}^{j+1} -2T_i^{j+1} +T_{i-1}^{j+1} \right)
% or
% {-\alpha T_{i+1}^{j+1} +\left(1+2\alpha \right)T}_i^{j+1} -\alpha T_{i-1}^{j+1} \approx T_i^j
% Now, for each time step, we write this equation at all the mesh points, solve the resulting linear system for all the temperatures at the new time step, and then repeat for the next step.
% The boundary condition at the origin can be addressed as follows:
% \begin{array}{l}
% \frac{\partial T}{\partial x}\left(x=0\right)\approx \frac{T_2^{j+1} -T_0^{j+1} }{2h}\\
% -k\frac{T_2^{j+1} -T_0^{j+1} }{2h}=\mathrm{q"}\\
% T_0^{j+1} =T_2^{j+1} +\frac{2\mathrm{hq"}}{k}
% \end{array}
% We can also rewrite our PDE approximation at the surface with
% \begin{array}{l}
% {-\alpha T_2^{j+1} +\left(1+2\alpha \right)T}_1^{j+1} -\alpha T_0^{j+1} \approx T_1^j \\
% {-2\alpha T_2^{j+1} +\left(1+2\alpha \right)T}_1^{j+1} \approx T_1^j +\frac{2\alpha \mathrm{hq"}}{k}
% \end{array}
% Now we have what we need to implement this algorithm as follows:

clear all
q=1e5;
k=100; rho=1000; cp=500;
diffus=k/rho/cp;
length=1; 
ndiv=10; %number of spatial divisions
h=length/ndiv; %mesh spacing
endtime=300;
ntimes=150;
dt=endtime/ntimes;
alpha=diffus*dt/h^2
tboundary=0; %boundary temperature (C)
time=zeros(ntimes+1);
temperature=zeros(ndiv+1);
tsurf=zeros(ntimes+1);
told=zeros(ndiv+1);
a=(1+2*alpha)*eye(ndiv+1)-alpha*diag(ones(ndiv,1)',1)-alpha*diag(ones(ndiv,1)',-1);
%adjust for boundary condition at x=L
a(ndiv+1,ndiv+1)=1; a(ndiv+1,ndiv)=0;
%adjust for boundary condition at x=0
a(1,2)=-2*alpha;
for i=1:ntimes
    b=temperature;
    b(1)=b(1)+2*alpha*h*q/k;
    b(ndiv+1)=0;
    time(i+1)=i*dt;
    temperature=a\b;
    tsurf(i+1)=temperature(1);
end
tanalytic=2*q/k*sqrt(time*k/rho/cp/pi);
plot(time,tsurf,time,tanalytic,'o')

% Here I've taken 1/10 as many time steps as before, but still obtained reasonable results.
% © James Blanchard, University of Wisconsin-Madison, blanchard@engr.wisc.edu