This lecture starts off with Example 15.4 from the notes, where we look to solve $$ \begin{gathered} u_t = u_{xx} \\ u(0, t) = 2, \qquad u(\pi, t) = 1 \\ u(x, 0) = 0. \end{gathered} $$
We show that the solution is given by $$ u(x, t) = U(x) + \sum_{n=1}^\infty B_n \sin(nx) \mathrm{e}^{-n^2 t} $$ where we have found the steady-state solution $$ U(x) = 2 - \frac{x}{\pi}, $$ as well as the coefficients $$ B_n = \frac{2}{n\pi}[(-1)^n - 2]. $$
We also illustrated the solution using the code:
Inhomogeneous heat equation
%% Plot the Fourier series for Example 15.4 % Written for MA20223 Vectors & PDEs 2019-20 clear % Clear all variables close all % Close all windows N = 100; % How many Fourier modes to include? R = @(n, t, x) 2/(n*pi)*((-1)^n - 2)*exp(-n^2*t)*sin(n*x); % Create a mesh of points between two limits x0 = pi; x = linspace(0, x0, 1000); % Create a mesh of points in time t = linspace(0, 5, 200); figure(1); % Open the figure plot(x, 2 - x/pi, 'b', 'LineWidth', 2); % Plot the base function ylim([-0.2,2]); % Set the y limits xlim([0, x0]); % Set the x limits xlabel('x'); ylabel('u(x,t)'); hold on for j = 1:length(t) tj = t(j); u = (2 - x/pi); for n = 1:N u = u + R(n, tj, x); end if j == 1 p = plot(x, u, 'r'); else set(p, 'YData', u); end drawnow title(['t = ', num2str(tj)]); pause(0.1); if j == 1 pause end end
The next thing we studied was PS8, Q1, which is the solution of the homogeneous Neumann problem for the heat equation.
Modification of PS8 Q1
%% Plot the Fourier series for a made-up modification of PS8 Q1. % Written for MA20223 Vectors & PDEs clear % Clear all variables close all % Close all windows N = 20; % How many Fourier modes to include? % Define an in-line function that takes in three inputs: R = @(n, t, x) 4*(-1)^n/n^2*cos(n*x)*exp(-n^2*t); % Create a mesh of points between two limits x0 = pi; x = linspace(0, x0, 1000); % Create a mesh of points in time t = linspace(0, 5, 200); figure(1); % Open the figure plot(x, x.^2, 'b', 'LineWidth', 2); % Plot the base function ylim([-0.2,pi^2]); % Set the y limits xlim([0, x0]); % Set the x limits xlabel('x'); ylabel('u(x,t)'); hold on for j = 1:length(t) tj = t(j); u = 1/2*(2*pi^2/3); for n = 1:N u = u + R(n, tj, x); end if j == 1 p = plot(x, x.^2, 'r'); else set(p, 'YData', u); end drawnow title(['t = ', num2str(tj)]); pause(0.1); if j == 1 pause end end