This lecture starts off with Example 15.4 from the notes, where we look to solve ut=uxxu(0,t)=2,u(π,t)=1u(x,0)=0.
We show that the solution is given by u(x,t)=U(x)+∞∑n=1Bnsin(nx)e−n2t where we have found the steady-state solution U(x)=2−xπ, as well as the coefficients Bn=2nπ[(−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