%% Numerical Solution of the Heat Equation
%
% In this section, we will use MATLAB to numerically solve the heat 
% equation (also known as the diffusion equation), a partial 
% differential equation that describes many physical processes 
% such as conductive heat flow or the diffusion of an impurity in a 
% motionless fluid.  You can picture the process of diffusion as a drop 
% of dye spreading in a glass of water.  (To a certain extent you could 
% also picture cream in a cup of coffee, but in that case the mixing is 
% generally complicated by the fluid motion caused by pouring the cream 
% into the coffee, and is further accelerated by stirring the coffee.) 
% The dye consists of a large number of individual particles, each of 
% which repeatedly bounces off the surrounding water molecules, 
% following an essentially random path.  There are so many dye 
% particles that their individual random motions form an essentially 
% deterministic overall pattern as the dye spreads evenly in all 
% directions (we ignore here the possible effect of gravity).  In a 
% similar way, you can imagine heat energy spreading through random 
% interactions of nearby particles.
%
% In a three-dimensional medium, the heat equation is
%
% $$\frac{\partial u}{\partial t}=k\left( \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2} + \frac{\partial^2u}{\partial z^2}\right).$$
%
% Here $u$ is a function of $t$, $x$, $y$, and $z$ that represents the 
% temperature, or concentration of impurity in the case of diffusion, 
% at time $t$ at position $(x, y, z)$ in the medium.  The constant $k$
% depends on the materials involved, and is called the thermal 
% conductivity in the case of heat flow, and the diffusion coefficient 
% in the case of diffusion.  To simplify matters, let us assume that 
% the medium is instead one-dimensional.  This could represent 
% diffusion in a thin water-filled tube or heat flow in a thin 
% insulated wire; let us think primarily of the case of heat flow.  
% Then the partial differential equation becomes
% 
% $$\frac{\partial u}{\partial t}=k\,\frac{\partial^2u}{\partial x^2},$$
% 
% where $u(x, t)$ is the temperature at time $t$ a distance $x$ along 
% the wire.

%% A Finite-Difference Solution
%
% To solve this partial differential equation, we need both initial 
% conditions of the form $u(x, 0) = f(x)$, where $f(x)$ gives the 
% temperature distribution in the wire at time 0, and boundary 
% conditions at the endpoints of the wire; call them $x = a$ and $x = b$.  
% We choose so-called Dirichlet boundary conditions $u(a, t) = T_a$
% and $u(b, t) = T_b$, which correspond to the temperature being held 
% steady at values $T_a$ and $T_b$ at the two endpoints.  Though an 
% exact solution is available in this scenario, let us instead 
% illustrate the numerical method of finite differences.
%
% To begin with, on the computer we can keep track of the 
% temperature $u$ only at a discrete set of times and a discrete set of 
% positions $x$.  Let the times be $0, \Delta t, 2\Delta t,\dots, N\Delta t$, 
% and let the positions be $a, a + \Delta x,\dots, a + J\Delta x = b$, 
% and let $u^n_j = u(a + j\Delta t, n\Delta t)$.  After rewriting the partial 
% differential equation in terms of finite-difference approximations to 
% the derivatives, we get
% 
% $$\frac{u^{n+1}_j-u^n_j}{\Delta t}=k\frac{u^n_{j+1}-2u^n_j+u^n_{j-1}}{\Delta
% x^2}.$$
% 
% (These are the simplest approximations we can use for the derivatives,
% and this method can be refined by using more accurate approximations, 
% especially for the $t$ derivative.)  Thus if, for a particular $n$, we 
% know the values of $u^n_j$ for all $j$, we can solve the equation 
% above to find for each $j$,
% 
% $$u^{n+1}_j=u^n_j+\frac{k\Delta t}{\Delta x^2}(u^n_{j+1}-2u^n_j+u^n_{j-1})
% =s(u^n_{j+1}+u^n_{j-1})+(1-2s)u^n_j,$$
%
% where $s = k\Delta t/(\Delta x)^2$.  In other words, this equation 
% tells us how to find the temperature distribution at time step $n+1$ 
% given the temperature distribution at time step $n$.  (At the 
% endpoints $j = 0$ and $j = J$, this equation refers to temperatures 
% outside the prescribed range for $x$, but at these points we will 
% ignore the equation above and apply the boundary conditions instead.)
% We can interpret this equation as saying that the temperature at a 
% given location at the next time step is a weighted average of its 
% temperature and the temperatures of its neighbors at the current time 
% step.  In other words, in time $\Delta t$, a given section of the wire
% of length $\Delta x$ transfers to each of its neighbors a portion $s$ 
% of its heat energy and keeps the remaining portion $1-2s$.
% Thus, our numerical implementation of the heat equation is a 
% discretized version of the microscopic description of diffusion we 
% gave initially -- that heat energy spreads due to random interactions 
% between nearby particles.
%
% The following M-file, which we have named |heat.m|, iterates the 
% procedure described above.

type heat

%%
% The function *|heat|* takes as inputs the value of $k$, vectors of $t$- and
% $x$-values, a vector *|init|* of initial values (which is assumed to have 
% the same length as *|x|*), and a vector *|bdry|* containing a pair of 
% boundary values.  Its output is a matrix of $u$-values.  Notice that 
% since indices of arrays in MATLAB must start at $1$, not $0$, we have 
% deviated slightly from our earlier notation by letting $n=1$ represent 
% the initial time and $j=1$ represent the left endpoint.  Notice also 
% that, in the first line following the *|for|* statement, we compute an 
% entire row of *|u|*, except for the first and last values, in one line; 
% each term is a vector of length $J-2$, with the index $j$ increased by $1$ 
% in the term *|u(n,3:J)|* and decreased by $1$ in the term *|u(n,1:J-2)|*.
%
% Let's use the M-file above to solve the one-dimensional heat equation 
% with $k = 2$ on the interval $-5 \leq x \leq 5$ from time $0$ to time $4$,
% using boundary temperatures $15$ and $25$, and an initial temperature 
% distribution of $15$ for $x < 0$ and $25$ for $x > 0$.  You can imagine 
% that two separate wires of length $5$ with different temperatures are 
% joined at time 0 at position $x = 0$, and each of their far ends 
% remains in an environment that holds it at its initial temperature.  
% We must choose values for $\Delta t$ and $\Delta x$; let's try 
% $\Delta t = 0.1$ and $\Delta x = 0.5$, so that there are $41$ values 
% of $t$ ranging from $0$ to $4$ and $21$ values of $x$ ranging from $-5$ to
% $5$.

tvals = linspace(0, 4, 41);
xvals = linspace(-5, 5, 21);
init = 20 + 5*sign(xvals);
uvals = heat(2, tvals, xvals, init, [15 25]);
surf(xvals, tvals, uvals)
xlabel x; ylabel t; zlabel u  

%%
% Here we used *|surf|* to show the entire solution $u(x, t)$.  The output 
% is clearly unrealistic; notice the scale on the $u$-axis!  The 
% numerical solution of partial differential equations is fraught with 
% dangers, and instability like that seen above is a common problem 
% with finite-difference schemes.  For many partial differential 
% equations, a finite-difference scheme will not work at all, but for 
% the heat equation and similar equations, it will work well with proper 
% choice of $\Delta t$ and $\Delta x$.  We might be inclined to think 
% that, since our choice of $\Delta x$ was larger, it should be reduced, 
% but in fact this would only make matters worse.  Ultimately the only 
% parameter in the iteration we're using is the constant $s$, and one 
% drawback of doing all the computations in an M-file as we did above 
% is that we do not automatically see the intermediate quantities it 
% computes.  In this case we can easily calculate that
% $s = 2(0.1)/(0.5)2 = 0.8$.
% Notice that this implies that the coefficient $1-2s$ of $u^n_j$ in 
% the iteration above is negative.  Thus, the ``weighted average'' we 
% described before in our interpretation of the iterative step is not a 
% true average; each section of wire is transferring more energy than 
% it has at each time step!
%
% The solution to the problem above is thus to reduce the time step 
% $\Delta t$; for instance, if we cut it in half, then $s = 0.4$, and 
% all coefficients in the iteration are positive:

tvals = linspace(0, 4, 81);
uvals = heat(2, tvals, xvals, init, [15 25]);
surf(xvals, tvals, uvals)
xlabel x; ylabel t; zlabel u  

%%
% This looks much better!  As time increases, the temperature 
% distribution seems to approach a linear function of $x$.  Indeed,
% $u(x, t) = 20 + x$ is the limiting ``steady state'' for this problem; 
% it satisfies the boundary conditions and it yields 0 on both sides of 
% the partial differential equation.
%
% Generally speaking, it is best to understand some of the theory of 
% partial differential equations before attempting a numerical solution 
% as we have done here.  However, for this particular case at least, 
% the simple rule of thumb of keeping the coefficients of the iteration 
% positive yields realistic results.  A theoretical examination of the 
% stability of this finite-difference scheme for the one-dimensional 
% heat equation shows that indeed any value of $s$ between $0$ and $0.5$
% will work, and suggests that the best value of $\Delta t$ to use for 
% a given $\Delta x$ is the one that makes $s = 0.25$.  Notice that,
% while we can get more accurate results in this case by reducing $\Delta x$,
% if we reduce it by a factor of $10$, we must reduce $\Delta t$ by a 
% factor of $100$ to compensate, making the computation take $1000$ times 
% as long and use $1000$ times the memory!

%% The Case of Variable Conductivity
% Earlier we mentioned that the problem we solved numerically could also
% be solved analytically.  The value of the numerical method is that it 
% can be applied to similar partial differential equations for which an 
% exact solution is not possible or at least not known.  For example, 
% consider the one-dimensional heat equation with a variable coefficient,
% representing an inhomogeneous material with varying thermal 
% conductivity $k(x)$,
% 
% $$\frac{\partial u}{\partial t}=\frac{\partial}{\partial x}\left( k(x)\frac{\partial u}{\partial x}\right)=k(x)\frac{\partial^2u}{\partial x^2}+k'(x)\frac{\partial u}{\partial x}.$$
%
% For the first derivatives on the right-hand side, we use a symmetric 
% finite-difference approximation, so that our discrete approximation 
% to the partial differential equations becomes
% 
% $$\frac{u^{n+1}_j-u^n_j}{\Delta t}=k_j\frac{u^n_{j+1}-2u^n_j+u^n_{j-1}}{\Delta x^2}+\frac{k_{j+1}-k_{j-1}}{2\Delta x}\frac{u_{j+1}-u_{j-1}}{2\Delta x},$$
%
% where $k_j=k(a + j\Delta x)$.  Then the time iteration for this method is
% 
% $$u^{n+1}_j=s_j(u^n_{j+1}+u^n_{j-1})+(1-2s_j)u^n_j+0.25(s_{j+1}-s_{j-1})(u^n_{j+1}-u^n_{j-1}),$$
%
% where $s_j = k_j\Delta t/(\Delta x)^2$.  In the following M-file, 
% which we called |heatvc.m|, we modify our previous M-file to 
% incorporate this iteration:

type heatvc

%%
% Notice that *|k|* is now assumed to be a vector with the same length 
% as *|x|*, and that as a result so is *|s|*.  This in turn requires that 
% we use vectorized multiplication in the main iteration, which we have 
% now split into three lines.
%
% Let's use this M-file to solve the one-dimensional variable-coefficient
% heat equation with the same boundary and initial conditions as before,
% using $k(x) = 1 + (x/5)^2$.  Since the maximum value of $k$ is $2$, we 
% can use the same values of $\Delta t$ and $\Delta x$ as before:

kvals = 1 + (xvals/5).^2;
uvals = heatvc(kvals, tvals, xvals, init, [15 25]);
surf(xvals, tvals, uvals)
xlabel x; ylabel t; zlabel u  

%%
% In this case, the limiting temperature distribution is not linear; it 
% has a steeper temperature gradient in the middle, where the thermal 
% conductivity is lower.  Again we could find the exact form of this 
% limiting distribution, $u(x, t) = 20(1 + (1/\pi)\arctan(x/5))$, by 
% setting the $t$ derivative to zero in the original equation and 
% solving the resulting ordinary differential equation.
%
% You can use the method of finite differences to solve the heat 
% equation in two or three space dimensions as well.  For this and 
% other partial differential equations with time and two space 
% dimensions, you can also use the PDE Toolbox, which implements the 
% more sophisticated ``finite-element method.''

%% A Simulink Solution
% We can also solve the heat equation using Simulink.  To do this, we 
% continue to approximate the $x$-derivatives with finite differences, 
% but think of the equation as a vector-valued ordinary differential 
% equation, with $t$ as the independent variable.  Simulink solves the 
% model using MATLAB's ODE solver, *|ode45|*.  To illustrate how to do this,
% let's take the example we started with, the case in which $k = 2$ on 
% the interval $-5\leq x\leq 5$ from time $0$ to time $4$, using boundary 
% temperatures $15$ and $25$, and an initial temperature distribution of
% $15$ for
% $x<0$ and $25$ for $x>0$.  We replace $u(x,t)$ for fixed $t$ by the vector
% *|u|* of values of $u(x,t)$, with, say, *|x = -5:5|*.  Here there are $11$ 
% values of $x$ at which we are sampling $u$, but since $u(x,t)$ is 
% predetermined at the endpoints, we can take *|u|* to be a nine-dimensional 
% vector, and we just tack on the values at the endpoints when we
% have finished.
% Since we're replacing $\partial^2/\partial x^2$ by its finite difference
% approximation and we've taken $\Delta x = 1$ for simplicity, our equation
% becomes the vector-valued ODE 
% 
% $$\frac{\partial{\mathbf u}}{\partial t}=k(A\mathbf u+c)$$
% 
% Here the right-hand side represents our approximation to 
% $k(\partial^2 u/\partial x^2)$.  The matrix $A$ is
%
% $$\begin{pmatrix}-2&1&\dots&0\\1&-2&\ddots&\vdots\\
% \vdots&\ddots&\ddots&1\\0&\dots&1&-2\end{pmatrix},$$
% 
% since we are replacing $\partial^2 u/\partial x^2$ at $(n, t)$ with
% $u(n-1,t) - 2u(n,t) + u(n+1, t)$.  We represent this matrix in MATLAB's
% notation by *|-2*eye(9) + diag(ones(8,1), 1) + diag(ones(8,1), -1)|*.
% The vector *|c|* comes from the boundary conditions, and has $15$ in its 
% first entry, $25$ in its last entry, and $0$'s in between. We represent it 
% in MATLAB's notation as *|[15; zeros(7,1); 25]|*.  The formula for *|c|*
% comes
% from the fact that *|u(1)|* represents $u(-4,t)$, and 
% $\partial^2 u/\partial x^2$ at this point is approximated by 
% 
% $$u(-5,t) - 2u(-4,t) + u(-3, t) = 15 - 2\mathbf{u(1) + u(2)}$$
%
% and similarly at the other endpoint.  Here's a Simulink model representing
% this equation:

open_system heateq

%%
% Note that we need to specify the initial conditions for $u$ as 
% Block Parameters for the Integrator block, and that in the Block 
% Parameters dialog box for the Gain block, we need to set the 
% multiplication type to ``Matrix''.  Since *|u(1)|* through *|u(4)|* represent
% the solution
% $u(x, t)$ at $x = -4$ through $-1$, and *|u(6)|* through *|u(9)|* represent
% $u(x, t)$ at $x = 1$ through 4, we take the initial value of $u$ to be 
% *|[15*ones(4,1); 20; 25*ones(4,1)]|*.  (We use $20$ as a compromise
% at $x = 0$, since
% this is right in the middle of the regions where $u$ is $15$ and $25$.)  
% The output from the model is displayed in the *Scope* block in the form of
% graphs of the various entries of $u$ as a function of $t$, but it's more 
% useful to save the output to the MATLAB Workspace and then plot it with 
% *|surf|*.   Incidentally, it helps to reset the $y$-axis limits
% on the *Scope* block to run from $15$ to $25$. To make this adjustment
% and run the model from $t=0$ to $t=4$, we execute the commands:

set_param('heateq/Scope', 'YMin', '15');
set_param('heateq/Scope', 'YMax', '25');
hout = sim('heateq', 'SaveTime','on','TimeSaveName','t', ...
   'SaveState','on','StateSaveName','u','StopTime','4');
tout = hout.find('t'); uout = hout.find('u');
simplot(tout, uout)
set(gcf, 'InvertHardcopy', 'off')
%%
% This saves the simulation times (from $0$ to $4$)
% as a vector *|tout|* and the computed values of $u$ as *|uout|*, a matrix
% with nine columns.  Each row of these 
% arrays corresponds to a single time step, and each column of *|uout|*
% corresponds to one value of $x$.  But remember that we have to add in 
% the values of $u$ at the endpoints as additional columns in $u$.  So we 
% plot the data as follows:

u = [15*ones(length(tout),1), uout, ...
    25*ones(length(tout),1)];
x = -5:5;
clf reset
set(gcf, 'Color', 'White')
surf(x, tout, u)
xlabel('x'), ylabel('t'), zlabel('u') 
title('Solution to heat equation in a rod')  

%%
% Notice how similar this is to the picture obtained before for constant
% conductivity $k = 2$.  We leave it to 
% the reader to modify the model for the case of variable conductivity.

%% Solution with pdepe
% MATLAB has a built-in solver *|pdepe|* for partial 
% differential equations in one space dimension (as well as time $t$).
% To find out more about it, read the online help on *|pdepe|*.  The 
% instructions for use of *|pdepe|* are quite explicit but somewhat 
% complicated.  The method it uses is somewhat similar to that used in 
% the Simulink solution above; i.e., it uses an ODE solver in $t$ and 
% finite differences in $x$.  The following M-file solves the second 
% problem above, the one with variable conductivity. Note the use of
% function handles and subfunctions:

type heateqex2

%% 
% Running it gives:

heateqex2  

%%
% Again the results are very similar to those obtained before.
