Solving the Heat Equation

Take a metal bar which is hot at one end (suppose it was dipped in boiling water for a few seconds). What happens to it next? The heat moves down the bar--the hot end gets cooler and the cold end gets warmer--until the entire bar is a uniform temperature. Finally, after a long time, the entire bar will cool down to the same temperature as its surroundings.

The exact mathematics for how this process takes place was first described by Joseph Fourier in 1822, and modified by many others in the 200 years since. Today, we know it as the heat equation, and it is a cornerstone of many engineering simulations today.

In this assignment, you will implement a very simple simulation of the heat equation. Fully understanding how to work with the heat equation and its variants could be several lifetimes of work!

Background

The simplest algorithm for solving the heat equation on a 2D grid works by determining how much each cell should heat up or cool down based on its neighbors. We then heat or cool each cell by a small amount, and repeat this many times.

Intuitively, if a cell's is cooler than its neighbors, it should warm up, because heat from its neighbors flows into it. Conversely, if a cell is warmer than its neighbors, it will cool down, because its heat flows into its neighbors.

This intuition is captured by our heating rule: if we are looking at cell $i,j$, then the amount it should heat up or cool down during the next timestep is:

$$ \bigg(T[i +1 ,j] - T[i,j]\bigg) + \bigg(T[i-1, j] - T[i,j]\bigg) + \bigg(T[i,j+1] - T[i,j]\bigg) + \bigg(T[i,j-1] - T[i,j]\bigg) $$

Note that I have grouped these terms: the first term is the difference in temperature between the cell's right neighbor and itself, the second term is the difference in temperature between the cell's left neighbor and itself, etc.

We can get rid of the parentheses and rename the terms in order to get the following equation:

$$ \Delta[i,j] = T[i +1 ,j] + T[i-1, j] + T[i,j+1] + T[i,j-1] - 4T[i,j] $$

We will use this as our rule for how to heat/cool cells.

This rule is incomplete, however: how should we compute heat change for cells on the edge of the square, where some of these neighbors might not exist? For example, if we are at LaTeX: (0,0) , the LaTeX: i-1 cell does not exist!

To solve this problem we'll add a rule that says that all cells that are on the edge of the grid have temperature zero. You can imagine this as submerging the grid into a bath made of ice water--the water keeps the material at the very edge at exactly 0C.

Finally, we may want to control how fast the simulation runs. In order to do this, we will use a "timestep", which simulates how much time passes every time we do one of these steps. The actual change in temperature in the cell should be the computed $LaTeX: \Delta$ times the timestep.

Assignment

Write a function heat_equation that asks the user for an input file. The program will read all of its input from that file and execute the heat equation simulation on the result.

The input file format will be similar to what we used in Minesweeper. However, there will be a few additional numbers at the top of the file. In order, they are:

The timestep to use for the simulation
The number of timesteps to run
The number of rows/columns in the grid
The starting grid for the simulation

The first two numbers are floats, the third is an int, and the last is a grid of floating point numbers. You may assume that the edge conditions (all cells on the edge of the grid have temperature 0.0) will be true in your input.

An example input file is given at the end of this assignment.

Your program should read the input file, store the timestep and number of steps, and then run the heat equation simulation. At each step, you should compute the $\Delta$ given by the heating rule, then multiply this by the timestep. Once all $LaTeX: \Delta$ values have been computed, add them to the grid to execute a timestep.

Repeat this until you have executed all the requested timesteps, then print the final grid to the console.

IMPORTANT: You should not modify the grid of temperatures until after you have computed all the $LaTeX: \Delta$ values!

Insight

Take the sample file from the end of the assignment and run it for many, many timesteps. What happens to the temperatures in the grid? Remembering that the boundary conditions we chose correspond to being submerged in an ice water bath, does this behavior make sense?

Submission

Submit a single file named Heat.py on Canvas. Your file needs to compile and run. It should also have a header with the following information (this goes in your source file, not in the program output):

# File: Heat.py
# Student: 
# Course: Intro to Programming
# 
# Date:
# Description of Program:
# Explanation of behavior of program in the long-term:

The description should be a short (1-3 sentence) description of what the program does. Do not describe how it's written!

Reference Program

Since this code is very hard to debug, I have provided a sample program for you to run. This program consists of compiled python bytecode. You can run it the same way you would run a normal python file (e.g. python3 sample.pyc).

The reference program will ask one extra question compared to your program: it will ask you which frame you want to view. As opposed to your program, where you run until the end of the simulation and then print the last grid, the reference program can print snapshots in the middle of the simulation. You can use this to debug your simulation.

Note that "0" is the grid at before any updates have been applied, i.e. the one you read from the file.

Example File

0.1
1000
10
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 9.0 9.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 9.0 9.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0