Question 1. The 2D diffusion equation u_t = D∇²u is discretized using the Finite Difference method as

u_i,jⁿ⁺¹ - uⁿ_i,j/DΔt = γ [(uⁿ⁺¹_i+1,j - 2u_i+jⁿ⁺¹ + u_i-1,jⁿ⁺¹)/Δx² +  (uⁿ⁺¹_i,j+1 - 2u_i,jⁿ⁺¹ + u_i,j-1ⁿ⁺¹)/Δy²] +(1-γ)[(uⁿ_i+1,j - 2u_i+jⁿ + u_i-1,jⁿ)/Δx² +  (uⁿ_i,j+1 - 2u_i,jⁿ + u_i,j-1ⁿ)/Δy²],

where γ is a parameter between 0 and 1 which determines the method of time discretization. Specifically, we have

γ = 0 : explicit method

γ = 1 : implicit method

γ = 1/2 : Crank-Nicholson method

Derive the stability condition for each of the above time discretization methods.

Question 2. The propagation of electromagnetic waves in 2D is governed by the wave equation

n²/c².∂²u/∂t² = ∂²u/∂x² + ∂²u/∂y²,

where n(x, y) is the refractive index of the medium and c = 3 × 10⁸ m/s is the speed of light in vacuum. We would like to solve the above equation using the Finite Difference Time Domain (FDTD) method. The computation domain is restricted to a rectangular region of size a × b and Radiating Boundary Conditions based on the one-way wave equation are applied to all four boundaries.

(a) Give a FD discretization of the wave equation for an interior node (i, j).

(b) Derive the FD equations for nodes on the left, right, top and bottom boundaries.

(c) Implement the FDTD method in a MATLAB program to solve the above wave equation. An outline of the program is given at the end of the assignment for your reference.

(d) Use your program to run the following simulations:

Simulation 1:

This is a test run to verify that your program works correctly. Define the computational domain to be the region 0 ≤ x ≤ 10μm, -5μm ≤ y ≤ 5μm, and set the index of the medium to be free space everywhere, n(x, y) = 1. Apply a point source located at (2.5μm, 0). The source emits a Gaussian pulse modulating a carrier signal given by:

u(x_s, y_s, t) = exp[- ((t - T_o)/(w/2))²]sin(ωt) ,

where ω = 2Πc/λ is the frequency, λ = 1.0μm is the wavelength, w = 8.0fs is the pulse width, and T₀ = 4.0fs is the time offset of the pulse centre. Use grid sizes Δx = Δy = 0.05μm and set the time step to the maximum allowable by the CFL stability condition.

Initialize the fields at the first two time steps to 0 (i.e., set u⁰_i,j = u¹_i,j = 0). Run the simulation for 175 time steps and provide the following results:

(i) A 3D plot or a contour plot of the field distribution u(x, y) at time step n = 175. (Use mesh or surf commands for 3D plot and contour for contour plot)

(ii) A plot of the field u versus x along the line y = 0 at time step n = 175. Print out and hand in a copy of your MATLAB program.

Simulation 2:

In this simulation we would like to study the diffraction of a plane wave by a single slit. We consider a plate with a slit of width a = 5μm illuminated by a plane wave with wavelength λ = 1.0μm, as shown in the figure below. The plate is 0.5μm thick and is assumed to be made of a material of very high refractive index (n = 10) so that it is strongly reflective. In the simulation, set the computational domain to be 0 ≤ x ≤ 25μm, -15μm ≤ y ≤ 15μm, with grid sizes Δx = Δy = 0.05μm. Use a line source located at xs = 1.0μm to generate a plane wave with unit amplitude and wavelength λ = 1.0μm. The field at a node on the line source is given by u(xs , y, t) = sin(ωt) where -14.5μm ≤ y ≤ 14.5μm (so that the line source does not touch the top and bottom boundaries). Set the time step to the maximum allowable by the CFL stability condition, run the simulation and provide the following results:

(i) A 3D plot or a contour plot of the field distribution u(x, y) at t = 90fs.

(ii) Suppose a "screen" is placed at location xo = 20μm to image the diffraction pattern of the field intensity. The field distribution on the screen can be expressed as

u(x_o, y, t) = E(x_o, y) cos(ωt),

where E(x_o, y) is the envelope of the field. Give a plot of the intensity diffraction pattern, I(y) = E2(x_o, y). (The field envelope can be obtained by capturing the maximum field at each point on the screen over a few periods of oscillation of the wave).

(iii) From Fraunhofer's theory of diffraction, the intensity distribution of light transmitted through a single slit is given by

I(y) = I₀ [sin(β)/β]²

where β = (Πa / λ) sin θ. The angle θ is given by tan θ = y/L , where L = 10μm is the distance from the plate to the screen. Choosing an appropriate value for I_o, plot the above expression on the same plot obtained in (ii) and compare the numerical simulation to the theoretical intensity distribution.

Outline of the 2D Finite Difference Time Domain Program

% Create and initialize Ny x Nx matrices to store fields at time levels
% n - 1, n, and n + 1:

En_1 = zeros(Ny,Nx); % time level n - 1 En0 = zeros(Ny,Nx); % time level n
En1 = zeros(Ny,Nx); % time level n + 1

% March solution from time step n = 0 to n = nstop: for n = 1 : nstop,
% Compute interior nodes: for jj = 2 : Ny-1,
for ii = 2 : Nx-1, compute En1(jj,ii)
end
end

Compute value of E at source point (xs, ys) at time step n + 1 Compute boundary nodes at time step n + 1
% Store fields at the two most recent time levels:
Un_1 = Un0; Un0 = Un1;

end