Write a function to solve the two-dimensional in Matlab,unsteady heat conduction equation with no internal heat generationon a square domain. The length of each side of the square is 1 m.The function will have the following 4 inputs:

npts                    number of grid points in each coordinate direction

nt                        number of time steps to take

dt                        size of time step (s)

alpha                  thermal diffusivity (m²/s)

Use the following initial and boundary conditions:

Initialize the body to T = 0 ^oC. At time t = 0,boundary conditions are imposed such that the temperatures on thesides of the square vary linearly. Here are the cornertemperatures:

bottom left corner:         T=   50 ^oC

bottom right corner:      T = 500^oC

top rightcorner:            T = 350 ^oC

top left corner:             T = 150 ^oC

Write the function definition statement. If you type docfunction at the MATLAB command prompt it will give youdetails on the proper syntax. The syntax is general, showing how touse both inputs and outputs. We will have no outputs, only the fourinputs noted above.

Make the name of your function the same as the name of yourfile. For example, if you call your function heat,then your filename should be heat.m.

Create a variable name for each of the four corner boundaryconditions and the initial condition, and set their values.

Noting that the temperature on each edge varies linearly,calculate the temperature at every grid point on the outerboundaries. You’ll need the corner temperatures to do thiscalculation.

Create an array to contain the “old” temperature field at timelevel p, and call it told. A convenient way to do this is to usethe ones function (look at doc ones), whichcreates an array of all ones of specified dimensions. An easy wayto set the field to the initial condition is to multiply this arrayby the initial temperature.

If there are any calculated variables that you will use multipletimes, calculate them. For example, it’s handy to have a variablenpm = npts – 1, where npts is the input number of points in eachdirection.

Create an array to contain the “current” solution at time levelp+1, and call it tnew. An easy way to do initialize it in MATLAB isto simply use

tnew = told;

Create a time loop going from 1 to nt, where nt is the inputnumber of time steps.

Inside the time loop, create two other nested loops, one for the“m” indices (x direction) and one for the “n” indices (ydirection). Since we don’t have to re-calculate any values onboundaries, the loops will go from 2 to npm.

For each pair of indices within the loops, calculate thesolution at time level p+1 (tnew) as a function of temperatures attime level p (told) using the finite-difference form of the 2-Dheat equation.

We are going to do a simple animation of the unsteady results asthey evolve. Use the contourf function to create the evolvingcontour plot. (doc contourf) One of the arguments is the number ofcontours, and 20 seems to work well for this problem. We want tore-draw the contour plot once per time step, so contourf should becalled after the end of the nested spatial loops, but before theend of the time loop. It’s not necessary to use a plot handle as wedid in the example in class. Simply call contour each time throughthe loop.

MATLAB tries to be efficient, so when it sees the contourffunction, it will buffer the results rather than plotting themright away. We don’t want it to buffer since we want to create areal-time animation. To prevent buffering, include a drawnowstatement on a separate line immediately after the call tocontourf. (This is similar to the example we did in class.)

Now that the current time step is complete, we have to copy thenew solution to the old one so we can go on to the next iteration,told = tnew. Make this the last line before the end of the timeloop.

Add the following to the end of the time loop. The first lineadds a “colorbar,” which is a scale showing the correspondencebetween temperatures and colors. The second line adds a label tothe colorbar. Finally, the third line gets rid of the plot axes,since we don’t need them here.

hc = colorbar;

hc.Label.String = ‘Temperature, deg. C’;

axis off

When you’re got your program ready to run, use the followinginput values:

npts = 40

nt = 1000

dt = 0.1

alpha = 0.001

To run the program, simply type

heat(40, 1000, 0.1, 0.001)