Live Scripts for Online Teaching: Solving a Heat Equation Example

Learn how to use a Live Script to teach a comprehensive story about heat diffusion and the transient solution of the Heat Equation in 1-dim using Fourier Analysis:

• The Story: Heat Diffusion
  • The transient problem
  • The great Fourier’s ideas
  • Thermal diffusivity of different material
  • The Physics where The Heat Equation come from
  • Visualization decaying of the heat structure
  • Building Symbolic heat structure
  • Separation of variables with BC
  • Fourier Analysis of IC
  • Embedding Fourier coefficients
  • Develop a flexible local function to build the heat solution
  • Experimenting with different materials

  Published: 18 Feb 2021

  PART 1: Heat Diffusion (the Story)

  1. The Transient Problem [01:32]
     1. Hi, today I’m going to demonstrate how to teach a lecture about solving the Heat Equation using Live Script proficiently. In a Live Script you can describe a problem with some text, implement a solution with the code and then run to immediately get the output. Text, code and output altogether will create a compelling story. We learn from stories. If you need a quick review of Live Scripts features, please watch the video linked in the transcript.
     2. Now Let’s go deep into our lecture. What’s the problem? It’s heat diffusion along a thin 1-dim rod. An image can help visualize the problem immediately: given the initial temperature across the rod at time zero and holding boundaries at the same fixed temperature, we want to compute how the temperature will change due to the heat conduction across the rod and over time
     3. By getting deeper into the structure of the heat solution, we’ll be able to discover that transient state will exponentially decay toward a stationary equilibrium. The time decay will strongly depend on thermal behavior of material. Diffusivity will play a key role: this property describes the attitude of a material to transfer heat without heating up too quickly. All sections can be organized into a summary table that allows to scroll up and down or quickly jump to next section.
     1. To provide more background, we recall the key Fourier’s findings. Fourier analysis is the awesome idea to represent every function, even with jumps, as sum of harmonics. We’ll see how to compute Fourier coefficients of the initial temperature of the rod. Note that any equation can be easily added from Insert Tab and click on Equation to open the Equation Editor, where you can pick many symbols: sum, derivative, or integral.
     2. Fourier’s Law is the proportionality of heat flow rate and temperature gradient, the minus confirming the common experience that heat always flows from hot to cold irreversibly. The Heat Equation is an amazing concentration of Fourier’s Law, Thermodynamics and conservation of energy, where diffusivity pops up as the key parameter, this is the main equation we want to solve together.
     3. Finally, we have two important ratios to highlight: First, the relationship between D and k. Please don’t be confused: conductivity k tells us only about heat transfer rate, while diffusivity D considers also the temperature change occurred during the heat transfer. So larger diffusivity requires also small heat capacity. The other remarkable ratio is the time decay of the transient. So when changing the material, remember: the larger the diffusivity, the smaller the time decay. It means that the transient will decay faster.
     1. Time has come for an Exercise where we want to explore real numerical values of diffusivity. So using this external page, I found the values for different materials and I put them into a table and saved it into mat file. Load the mat file containing this table called materials. We can create a variable to map the material name directly to its corresponding diffusivity value. For example, we can take diffusivity for iron, ok 23 (mm^2/s). To help choosing the material name, we can open the table, copy the name column. Then go to Insert Tab, click Control and pick a Drop Down, paste the list of materials. Now everybody can scroll this list and pick any material without typing the name, for instance iron, steel, or copper or pyrolytic graphite. So, we have a quick way to determine the diffusivity for any material in this list. We’d better save this information for future reusage.
     1. Live script can also help understand the physics behind the heat equation. Conservation of energy is the starting point. Look at units because they can tell us What energy balance is about: time change of heat energy is due to incoming heat flux. And this balance is in Watt, or Watt per volume in the differential form. Thermodynamics and Fourier’s Law allow to express energy and flux in terms of temperature. By embedding these laws into energy balance, the heat equation follows immediately. Now we can see where Diffusivity comes from.
     2. As an exercise. You can try to derive again the heat equation using the Equation Editor. In another exercise, we can learn how to use symbolic variables and assign units that are visualized in the output too. This way, you can verify the unit of diffusivity ratio. You can also define partial derivative equations. And check if they are consistent, that is the 2 members have the same units.

     Part II: Transient Heat Structure

     1. We’ve completed the first Part about the Story and some Key Facts. Now we are ready to start the second Part with more code to visualize the transient structure, build the heat solution, finally simulate heat diffusion in different materials. Let’s go. To solve the Transient problem, we focus on homogeneous boundary conditions, it means that boundaries are held constant at 0 Kelvin degrees. We’ll freeze these BC into a general heat structure. The solution structure is a series with space and time well “SEPARATED”, sine and exponential being eigenvectors of differential operators; boundary conditions are satisfied due the value of lambda; the negative exponential decaying to zero, the greater the Diffusivity the faster the decay. So, we’re going now to visualize this solution.
     2. We start writing the separated product sin(x)* exp(-t), this is actually a function handle on x, t, n. Lambda and b are new function handles as well depending only on n. For now, we can guess any expression for b. Later we’ll replace b with some Fourier coeff. Now, fix n = 1, what happens is we get a new function handle depending just on x and t, so this is a surface we can visualize pretty easily. With a simple for loop we can overlap multiple terms for different values of n. This is useful just to familiarize a bit with the single terms. We can copy and paste the code because I’d like to feel the shape of their sum. So we just need to accumulate the contribution from every term the get the final heat structure. To generalize with N terms, we could Insert a Slider control to change the number of terms we want to sum. When you click, you get the heat surface and the drawnow function will create this nice animation.
     1. We can also use a symbolic approach to build the heat structure. Copy the previous function handles and let’s see a trick to convert them into symbolic functions: just move the argument list from the right to the left. We must define the required symbolic variables, and the symbolic functions too. Let me redefine b as a variable, not a function, and let’s make assumptions that can simplify to compute Fourier coefficients later: n to be integer, L and D to be positive. When you run, symbolic output is nice. Note that symbolic objects can be either SYMBOLIC variables or SYMBOLIC FUNCTIONS. Now To reuse this code, we’d better put it into a local function at the end of the Live Script. You just type function, give a name, define heat structure as the output. So now we have a quick way to build the symbolic heat structure with the name you like just by calling this local function anytime we need.
     2. We cannot visualize the symbolic heat object until we fix some values. Now, substituting a symbolic variable with any value would be easy. But what if L is already used to store the numerical value? The previous syntax cannot work because we must redefine L as a Symbol, so now we are also able to substitute the symbol L with the numerical L. We can embed multiple parameters like D and b_n with a single command. But this has remained a function still depending on n. To fix also n = 3, we don’t need subs operator, simply place 3 as function input, and finally we get the surface depending just on space and time. And at last we could visualize it. But where is the heat? The heat is the sum of many of these terms. So to get the heat, we could compute the sum of the first 3 terms with n = 1, 2, and 3, we can reuse the same code to visualize also the sum, maybe increasing the time duration, so we’ve learnt how to embed also numerical parameters into the general heat structure.

     Part III: Solving Heat Problem

     1. Now we want to understand the why and where the heat structure come from. We’re going to verify analytical methods to show that Not only is the heat structure A solution, but also that it’s THE solution of the problem. In this small exercise we verify that heat structure satisfies the Heat Equation. MATLAB will compute the partial derivatives for us. When we run, the output shows that the heat structure is indeed, a solution. Suppose we have defined the heat problem, but we want to look for a solution. The first step will be separation of variables considering only the boundary conditions; the second step will be Fourier analysis to include also the initial condition.
     2. Separation of variables is the idea that solution is product of two separated functions, each one depending on a single variable. We substitute this product both in the equation and the boundary conditions. And Live Script will show the output immediately. We can further separate the equation, putting all the terms depending on time on the left and those depending on space on the right. The only way for the two sides to be equal for any x and t is that they are constant, a negative constant to avoid trivial solutions, so let me write it as negative squared lambda.
     3. We get two separated ODEs, first order for time, second order for space, and that’s where exponential and sine show up! We separate the two sides here and equate to minus squared lambda. Then we solve the ODEs separately using the left boundary condition, and suitable constants. Now it’s the turn of the right boundary condition that will force lambda to depend on the rod length. MATLAB can determine lambda spectrum to be any integer multiple of pi over L. So, we substitute lambda, rename all the constants as b_n, and taking the final product, we have demonstrated the heat solution must have this structure. But what about bn? None of these will satisfy the initial condition yet. So any idea about how to choose bn?
     1. At time 0 The heat structure boils down to a series of sine, but what to do to make it equal to the initial temperature? Bn are still free degrees of freedom. At the beginning we recalled Fourier’s Analysis to represent any function as sum of harmonics. But What is there behind? the magic idea is orthogonality. Very simply, just multiply each member by a sin and then take integral. Now Let’s compute this integral which is the inner product of 2 harmonics. The output shows a zero when indexes are different. So, harmonics are orthogonal! Bingo! To satisfy the initial condition we have this beautiful formula by taking the integral inner product and making a projection along the harmonics. Let’s do some examples.
     2. Let’s start with simple function, let’s say phi(x) = x and repeat the magic Fourier’s formula, phi(x) times the harmonic, then the integral from 0 to L, some factor to normalize and here are the coefficients. To visualize their behavior a little bit, we can first evaluate some values and convert them into double, so we create a bar plot to see they are decreasing, with an alternating sign. We could also plot the whole Fourier Series with the first N harmonics and taking the symbolic sum from 1 to N. So, with this for loop, we increment the number of harmonics we’re going to sum and then overlap on the original phi to feel how well they can approximate phi closer and closer. We could add a button to run, repeat the loop and quickly experiment with any other initial function. What happens if the initial condition is a square wave, with a jump in the middle? More interestingly, Fourier series can approximate even the jump. As last example, What about a saw tooth? Fourier series performs very well even in this case, no matter how many jumps you have. Amazing, isn’t it?
     1. Now, we review all the steps together to compute the final solution. First, we compute Fourier coefficients of the initial condition. Then, define the heat structure by reusing the local function we defined before. The structure satisfies the boundary conditions. Now, embed Fourier coefficients to satisfy the initial condition too, … we embed length and diffusivity as well. Finally, to get the heat solution, we take the sum of terms from 1 to N.
     2. Of course, we can change the length or any other parameter. For example, we could take an initial condition jumping from 10 to 40 degrees. When you run all computation will be updated. We could also change material diffusivity or increase the number of terms. To visualize the final surface, we also fix a time duration, time decay could be an option. Here is the heat solution. If we change initial condition, heat surface is responsive. To create an animation of the surface, we’d better use a numerical approach where we define a meshgrid and recompute the numerical heat matrix at each time step.

     Part IV: Final Simulation

     1. To simulate efficiently we’re going to create a more flexible function to manage all problem input, visualize the solution, or enable animation for a certain time, 30 minutes or 2 hours, depending on material. Going back to the previous section, we copy the 4 steps solving the problem and scroll down to a new local function where to paste them in a more compact and reusable way. So here we have a good synthesis of all we have learnt to solve the heat equation. Heat solution is part of the output arguments. Remarkably, the function can compute the material diffusivity by reusing the map container built at the beginning and provide suitable unit correction. It also takes care to convert given duration into seconds. Finally, it can enable animation checking last input.
     2. Testing a function is fundamental: if we run with 2 inputs, we get into troubles, why? We need 4 inputs to compute the heat solution. So, we can correct quickly. So now it’s ok. What about 4 inputs? It does work, but coefficients are a bit weird. We’d better add a breakpoint and step into the function. We execute step by step, continue to next breakpoint until we discover that n was not assumed to be an integer. So, we stop debugging and make the correction. By the way, “Oranges” also notifies two lines where semicolons are missing. Now test again another syntax with graphic, … check b, and test animation. Now the function is ready for simulation.
     1. Our flexible interface can help a lot, but still needs many inputs. Changing values of each of them may be time-consuming but Live controls can speed up this task. We already created such a Control at the very beginning to select materials. We used there a Drop Down to collect the whole list of materials. Note that when you click on a control, it can run the whole section but, in this case, I do Nothing. For example, let’s pick iron. There are different types of controls, for the rod length we can use an Edit Box to enter a double value. Let’s say 1 meter. For the initial condition, we just write the expression, while a Slider is perfect to choose the number of Fourier terms. We add another Edit box for the duration … and a Drop Down to select the time units.
     2. We can also add a Button to start simulation when we’re ready. We see iron has a small diffusivity, and reaches equilibrium in about 3 hrs. So, simulation says to change time scale from minutes to hours to catch the transient. We can try an initial condition with a jump and increase the number of terms. Or pick another material like copper with larger diffusivity about 100, so transient dies down much faster in about 30 minutes. Finally, Pyrolytic graphite is very important material in many applications because its diffusivity is greater than 1000 and reaches equilibrium in 3 minutes. Reducing the length, time scale becomes even seconds. To conclude, we can export any figure or the whole Live Script in different formats. For example, we can see the final Lecture about Heat Equation exported into a pdf file. So, I hope you enjoyed, and you can now use Live Script to diffuse your ideas. Thanks for watching