%% 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.