## The aim

Imagine we have a square coaster upon which we place a hot mug of tea. We wish to model the distribution of temperature across the coaster over time. For the sake of simplicity we will model the coaster only in two dimensions and we take the initial temperature across the surface to be except for a circle where the rim of the bottom of the mug touches the coaster, at which the temperature is .

In this post I’m going to show how we can model the heat equation succinctly in F#. I’m going to consider the two-dimensional case and approximate the solution at discrete spatial mesh points and at discrete time periods.

We will also plot the results by mapping the temperature onto the brightness (i.e. a heat or intensity map).

## The mathematics

In two dimensions the heat equation – taking the size of the coaster to be 100mm square – is given by:

where represents the temperature at time and at coordinates .

We need to apply boundary conditions at the edges of the coaster. We will assume for simpliciy that the temperature along the edges of the coaster remains constant, that is:

We also need to set our initial conditions:

To model this in F# we are going to represent the surface of the coaster using a 100×100 matrix (the matrix class is included in the F# powerpack).

Using the Euler method we can convert our continuous differential equation into a discrete difference equation:

For some constant which represents the thermal conductivity of the surface. Note that here is a natural number representing discrete time values.

## Show me the code!

The F# code runs very close to the mathematics so it should be self-documenting (although I’ve added some comments for readability). Plotting the results is relatively straightforward: we normalize the temperatures and represent them as shades of grey, white being hottest and black being coolest.

#light open Microsoft.FSharp.Math open System.Drawing open System.Drawing.Imaging open System.Windows.Forms open Microsoft.FSharp.Collections open System.Linq // Flattens a 2D array into a sequence let array2D_to_seq arr = seq {for i in 0..Array2D.length1 arr - 1 do for j in 0..Array2D.length2 arr - 1 do yield arr.[i,j]} // Find maximum value in a matrix let max_value_in_matrix m = m |> Matrix.toArray2D |> array2D_to_seq |> PSeq.max // Normalizes a matrix so its maximum value is 1 let normalize_matrix m = m * (1.0/(max_value_in_matrix m)) let mug_diameter = 50.0 //mm let coaster_length = 100.0 //mm let tolerance = 5.0 //we're operating on discrete space so the rim of the mug needs to have some thickness let num_steps = 1000 //number of iterations to be modelled // Number of rows and columns in the matrix let rows = (int)coaster_length let cols = (int)coaster_length // Equation for a circle let circle r x y = (x-coaster_length/2.0)**2.0 + (y-coaster_length/2.0)**2.0 - (mug_diameter/2.0)**2.0 // Inital conditions function let initialValues (x:int) (y:int) = match x,y with | (x,y) when circle (mug_diameter/2.0) (float(x)) (float(y)) >= 0.0 && circle (mug_diameter/2.0) (float(x)) (float(y)) <= tolerance**2.0 -> 80.0 |_ -> 20.0 // Create matrix representing initial conditions let initialConditions = Matrix.init rows cols initialValues |> normalize_matrix let c = 0.6 //Thermal conductivity let delta_t = ((1.0) / 2.0*c)/2.0 //Time interval // Our difference equation let rec temp_at x y (o:float) (l:float) (r:float) (t:float) (b:float) = o + c * delta_t * (r+l+4.0*o+t+b) // Mapping matrix u(t) to u(t+1) let newMatrix (m:matrix) = m |> Matrix.mapi(fun i j temp -> match (i,j) with | (i,j) when i = 0 || j = 0 || i = rows-1 || j = cols-1 -> 0.0 //Boundary conditions |_ -> temp_at i j (m.[i,j]) (m.[i-1,j]) (m.[i+1,j]) (m.[i,j+1]) (m.[i,j-1])) // Recursive function to determine the temperatures at time t let rec heatmap_at t = match t with | 0 -> initialConditions |_ -> heatmap_at (t-1) |> newMatrix let format = Imaging.PixelFormat.Format24bppRgb let toBitmap (arr:Color[,]) = // Create the bitmap let image = new Bitmap(arr.GetLength(0),arr.GetLength(1),Imaging.PixelFormat.Format24bppRgb) for i=0 to image.Width-1 do for j=0 to image.Height-1 do image.SetPixel(i, j, (arr.[i,j])) done done image let intensityMap intensity = Color.FromArgb((int (intensity * 255.0)),(int (intensity * 255.0)),(int (intensity * 255.0))) let intensities = heatmap_at num_steps |> normalize_matrix |> Matrix.toArray2D |> Array2D.map intensityMap let heatBitmap = intensities |> toBitmap let form = new Form( Text = "F# Heat Map", Size = heatBitmap.Size) let pic_box = new PictureBox( BorderStyle = BorderStyle.Fixed3D, Image = heatBitmap, Size = heatBitmap.Size, Dock = DockStyle.Fill, SizeMode = PictureBoxSizeMode.StretchImage) form.Controls.Add( pic_box ) #if INTERACTIVE form.Show() #else Application.Run(form) #endif

## Let’s take it for a spin!

Here I have taken snapshots at discrete times 0, 50, 100, …, 1000 with .

Quite an impressive simulation for just 100 lines of code – including comments and white space!

## References

Parallel Numerical Solution of 2-D Heat Equation, Verena Horak and Peter Gruber.

The Heat Equation, Wikipedia

Euler method, Wikipedia

Nice example! Like you, I’m finding one of the things I like best about F# is the compactness of the solutions it allows. Think about it, you’ve not only modeled the physics in 100 lines of code, you’ve implemented the display!

Pingback: Rick Minerich's Development Wonderland : F# Discoveries This Week 05/02/2010

Hi,

really enjoyed this – I wish you could find more of this stuff.

Hi James…

Very nice example. It really goes to show how succinctly you can represent models in F#.

One thing you might want to think about – have you had a look at the “Units of measure” feature of F#. You can then replace your comments,

let mug_diameter = 50.0 //mm

with,

let mug_diameter = 50.0

The advantage? Well, F# will now check at compile time that any calculations you perform have the correct units. The code will run the same, but we have an added layer of protection from bugs.

Thanks Andy, good suggestion. Units of measure are something I haven’t had a proper chance to look at yet but I will definitely try and take a look at them for a future example.

Watch this space. 🙂

Pingback: The Morning Brew - Chris Alcock » The Morning Brew #592

Neat example, but this would be only around 10 lines of code in something like MATLAB…

Hi James. Thanks for your F# blog posts. I’m not bored wrkn on RPP. I wish I had code to show you. F# is a help; I’m learning the math as I go.

An energy-driven approach to linkage unfolding – Cantarella, Demaine, Iben… present a new alg for unfolding planar polygonal linkages without self-intersection based on following the gradient flow of a” repulsive” energy function.

RPP – Refolding Planar Polygons – Iben, O’Brien, Demaine – Discrete & Computational Geometry, 2009 – Springer

… paper describes an alg for generating a guaranteed intersection-free

interpolation sequence between any pair of compatible polygons. … builds on prior results from linkage unfolding

Thanks. I had a brief skim of the paper – looks very interesting. I look forward to seeing an F# implementation!