## Transcribed Text

QUESTION 4
[40 MARKS TOTAL]
Background
In this question, you will be solving the 2D heat equation defined in Eqn 3 in a rectangular domain of
x 1 and -0.5 s y 505 The domain has an initial temperature of zero and the boundary
conditions are given by
T(x.y)=exp(1-x-y2)sin(x+y²)
Eqn. 4
You are asked to discretize the rectangular domain into Nx grid points in the x-direction and Ny grid
points in the y-direction. Note that the first and last points in both the x- and y-directions are
boundary conditions, which means that there will be a total of (Nx-2)x (Ny 2) unknowns
Q4a
in your PDF for Question 1. outline the algorithm needed to solve Eqn 3 using the BTCS scheme from
time = o to time t., with a constant timestep, At (assume that to NtAt). Include the major steps
required to set the problem up (don't include minor steps like setting grid spacing etc.) This
algorithm must be in the form of pseudo code or a flow chart. Example of a pseudo code for solving
an ODE using the Euler method is shown in the smippet below.
Set up parameters of the problem
Determine stepsize and number of steps
for each time step
Solve the Euler method
Store calculated values in vector
end
Q4b
Write a Matlab function that sets up the Laplacian matrix L and the boundary condition vector BC
Your function header must be written as
function BC) - Mat Set up (Nx. Ny. xg. yg, alpha, acf)
where
Nx and Ny are the number of points in the x- and y-directions, respectively
.
xg and Y8 are column vectors containing the grid points of , and y, respectively.
alpha is a, i.e. the thermal diffusion coefficient
BCF is function describing the boundary condition (see Eqn 4)
Lis the [(Nx-2)*(Ny-2) (Nx 2)*(Ny-2)] Laplacian matrix
BC is column vector on the RHS that accounts for the boundary conditions
NOTE: You can define an anonymous function that converts a 2D node number (j.k) = k) to a
global node number for the unknowns. Your function will look like
nn=@(j.k.Nx) {some function of k, Nx):
(See Workshops 25-27 for examples of what this looks like). When setting up the matrix problem
the function no (or at least a similar process) can be used to reliably determine the row and column
number of each entry in the Laplacian matrix and the row number in the BC matrix
Q4c
Modily Lab06_Q4.m to salve the 2D Heat equation using the BTCS scheme (i.e. Implement your
algorithm from Q4a) together with the function you have written in Q4b. Set a 0.1 and time step
At 1s, and solve the PDE from t 0 to =50s for (Na, Ny) (21, 11), [41, 21), (61, 31) and (81, 41).
Using the command 'tic' and 'toc', determine how much time is spent setting the matrix problem
up. i.e. the time spent assembling the Laplacian matrix versus how much time is spent undertaking
the time-stepping solution. (NOTE: The time taken for plotting should not be included as part of the
time taken for time-stepping).
Print out the times for each (Nx, Ny) grid and write . few sentences comparing the times and
explain what they mean.
In separate figures, plot contours of your salution Including the boundary values for each (Nx, Ny).
You will have to reshape your 10 vector of unknowns into a 20 matrix and embed that into a
slightly bigger 20 matrix and add boundary conditions.
Q4d
Repeat Q4c with timestep of At - 50: $ (Just do this inside the same script file (Lab06_Q4.m) with a
loop for the 2 different time steps).
Is the solution stable or unstable in this case? Should it be? Write a comment to the command
window that compares your solution calculated with one large time step to that obtained in 50
smaller time steps.
NOTE: the following MATLAB code will contour a 20 grid of data values stored in Oplot, where xg is a
1-D vector that specifies the grid x-locations (length=Nx) and Yg is a 1-D vector that specifies the grid
locations [length=Ny)
figure
clev-linspace(0.25,1.251: Set the desired contour levels

This material may consist of step-by-step explanations on how to solve a problem or examples of proper writing, including the use of citations, references, bibliographies, and formatting. This material is made available for the sole purpose of studying and learning - misuse is strictly forbidden.

% Applies BTCS algorithm for given parameters and returns:

% 1. the solution in matrix form (QnMax) with added boundary conditions

% 2. grid points on x and y axis (xg, yg)

% 3. total times needed to form Laplace matrix and to solve linear system

% (mat_time, sol_time)

function [QnMat, xg, yg, mat_time, sol_time] = BTCS(alpha, dt, t0, tmax, Nx, Ny)

xg = -1:2/(Nx-1):1;

yg = -0.5:1/(Ny-1):0.5;

BCf = @(x,y) exp(1-x^2-y^2)*sin(x^2+y^2);

nL = (Nx-2)*(Ny-2);

% Starting values are all zeros, outside boundaries

Qn = zeros(nL,1);

% matrix formation time and linear system solution time

mat_time = 0;

sol_time = 0;

for t=t0+dt:dt:tmax

% form the Laplace matrix and its RHS side from boundary conditions

tic;...