Question 7.4: Write a MATLAB program to solve ∂²T/∂x² + ∂²T/∂y² = 0 (7.4.2...

We have already worked out the interior node equation in Eq. (7.4.20). This equation will be valid for all of the nodes 2 ≤ i, j ≤ n − 1 except for the central node. We can write a short MATLAB program to handle the interior:

\frac{T_{k+1} − 2T_k + T_{k−1}}{Δ^2} + \frac{T_{k+n} − 2T_k + T_{k−n}}{Δ^2} = 0                                  (7.4.20)

This program, like several that follow, modifies the matrix A to account for a new part of the problem. Thus, it takes A as an input and returns it as an output. Note that we use nested for loops to handle the global counter k.

When we use the above program for the interior nodes, we are writing the wrong equation for the hot spot, so we need to fix this row in A. As we can see in Fig. 7.11, we need to use an odd number of grid points n in each direction if we want to have a node at the center. This node will be at the middle of the local index i

i =\frac{n + 1}{2}                             (7.4.23)

and the local index j,

j =\frac {n + 1}{2}                         (7.4.24)

The global index is

k = i + (j − 1)n                                                          (7.4.25)

so the center node is at

k_{center} = \frac{n + 1}{2}+ \left(\frac{n + 1}{2} − 1\right)n                           (7.4.26)

and we have

T_{k,center} = 1                                        (7.4.27)

We can implement these equations in a second MATLAB file:

This program computes the location of the mid-point, and then sets up a new row for A called insert_eqn that is all zeros except for the diagonal entry. We also modify the vector b to have the intensity of the hot spot.

For the top nodes, k = n(n − 1) + 2, ... , n^2 − 1, we have

T_k = 0                                (7.4.28)

which is identical to what we saw in our representative example. For the remaining boundaries, we need to make fictitious nodes as indicated in Fig. 7.11 For the bottom nodes k = 2, ... , n − 1, we need to use a fictitious node T_{k−n}

\frac{T_{k+n} − T_{k−n}}{2 Δx} = 0             (7.4.29)

which gives us

T_{k+n} = T_{k−n}                                    (7.4.30)

We then rewrite Eq. (7.4.20) as

2T_{k+n} + T_{k+1} − 4T_k + T_{k−1} = 0                                              (7.4.31)

For the left nodes, the fictitious node condition is very similar,

T_{k+1} = T_{k−1}                     (7.4.32)

where T_{k−1} is the fictitious node, which gives us

T_{k+n} + 2T_{k+1} − 4T_{k} + T_{k−n} = 0                                                                  (7.4.33)
For the right nodes, the fictitious node condition is the same as Eq. (7.4.32), except that T_{k+1} is now the fictitious node. As a result, we get

T_{k+n} + −4T_k + 2T_{k−1} + T_{k−n} = 0                                               (7.4.34)
All of the boundary nodes are implemented in the following program:

While it is possible to write the left nodes in terms of the global index k, it is simpler to just increment over the interior rows i = 1, ... , n−2 and then use k = ni+1. You could also use i = 2, ... , n − 1 and k = n(i − 1) + 1. Similar to the left nodes, it is easiest to count through the right nodes using two counters, i = 2, ... , n − 1 and then k = ni.
We now need to handle the corners. For the upper corners, the Dirichlet condition takes precedence and we have

T_{n(n−1)+1} = 0                                      (7.4.35)

and

T_{n^2} = 0                                (7.4.36)

For the lower corners, we have no-flux conditions. If you look at the direction of the normal vector to the lower left corner (k = 1) in Fig. 7.11, you can see that we need to consider both

\frac{∂T}{∂x} =\frac{ T_2 − T_1}{Δx}                  (7.4.37)

and

\frac{∂T}{∂y} = \frac{T_{n+1} − T_1}{Δy}                                           (7.4.38)

In the latter, we have used forward finite differences for convenience; a centered finite difference solution will give us the same result in the end. The no-flux condition is

n · ∇T = n_x\left(\frac{T_2 − T_1}{x}\right)+  n_y\left(\frac{T_{n+1} − T_1}{y}\right) = 0                    (7.4.39)

Since the normal vectors n_x = n_y for our problem, we simply get

T_{n+1} + T_2 − 2T_1 = 0                                           (7.4.40)

The same logic applies to the lower right corner (k = n), where we have

T_{2n} + T_{n−1} − 2T_n = 0                        (7.4.41)

The corners are thus implemented by again modifying A:

You may have noticed that the programs interior, edges, and corner only modify the matrix A, while the program hotspot modifies A and b. The reason we only modify A in the first three programs is because these equations are homogeneous. If we have inhomogeneous boundary conditions, such as a non-zero temperature at one of the boundaries, then we would need to modify both A and b there as well. We will see this change in Section 7.5.

Once we have all of the equations, we just need to use a technique from Chapter 2 to solve the resulting system. The program below solves the problem using a banded solver in MATLAB:

This is a very compact program – it sets up the space for A and b, fills them in using the programs that we already wrote, and then solves the problem. This modular approach to the program makes modifications very easy, since we have compartmentalized all of the boundary conditions and the governing PDE into their own subfunctions. One very important point to note is that the program takes in the number of nodes on an edge, n, so the matrix and vector are of size n^2.
The output of our main program is a vector of values of the temperature at each node. To make plots, we should unpack the vector in a manner similar to what we did for coupled BVPs in Section 6.6:

Figure 7.12 shows the solution to the problem. Visualizing these solutions is sometimes tricky, since they are functions of two variables. This particular projection was chosen so that you can easily see the Dirichlet boundary condition at the upper boundary and how the heat diffuses from the hot spot down to the sink on the upper boundary. On the other walls, the no-flux condition leads to a relatively boring result, where the temperature profile is relatively flat in the rest of the domain. We encourage you to play around with the solution to this problem, in particular to see how the location of the hot spot affects the solution. It is easy to move the hot spot around by just changing the variable k_hot in the program. You can also add additional hot spots relatively easily by adding additional subfunctions like hotspot that change other nodes from normal interior nodes to hotspots. For example, Fig. 7.13 shows two hot spots, one at half the intensity of the second. Clearly, we can construct some very interesting temperature profiles from this program!