Trinh @ Bath

Lecture 24: Computation of the heat equation III

Example 15.4: Solution of the inhomogeneous Dirichlet problem

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)en2t where we have found the steady-state solution U(x)=2xπ, as well as the coefficients Bn=2nπ[(1)n2].

We also illustrated the solution using the code:

Inhomogeneous heat equation

  1. %% Plot the Fourier series for Example 15.4
  2. % Written for MA20223 Vectors & PDEs 2019-20
  3.  
  4. clear % Clear all variables
  5. close all % Close all windows
  6. N = 100; % How many Fourier modes to include?
  7.  
  8. R = @(n, t, x) 2/(n*pi)*((-1)^n - 2)*exp(-n^2*t)*sin(n*x);
  9.  
  10. % Create a mesh of points between two limits
  11. x0 = pi; x = linspace(0, x0, 1000);
  12.  
  13. % Create a mesh of points in time
  14. t = linspace(0, 5, 200);
  15.  
  16. figure(1); % Open the figure
  17. plot(x, 2 - x/pi, 'b', 'LineWidth', 2); % Plot the base function
  18. ylim([-0.2,2]); % Set the y limits
  19. xlim([0, x0]); % Set the x limits
  20. xlabel('x'); ylabel('u(x,t)');
  21. hold on
  22.  
  23. for j = 1:length(t)
  24. tj = t(j);
  25. u = (2 - x/pi);
  26. for n = 1:N
  27. u = u + R(n, tj, x);
  28. end
  29. if j == 1
  30. p = plot(x, u, 'r');
  31. else
  32. set(p, 'YData', u);
  33. end
  34. drawnow
  35. title(['t = ', num2str(tj)]);
  36. pause(0.1);
  37. if j == 1
  38. pause
  39. end
  40. end

Problem set 8 Q1

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

  1. %% Plot the Fourier series for a made-up modification of PS8 Q1.
  2. % Written for MA20223 Vectors & PDEs
  3.  
  4. clear % Clear all variables
  5. close all % Close all windows
  6. N = 20; % How many Fourier modes to include?
  7.  
  8. % Define an in-line function that takes in three inputs:
  9. R = @(n, t, x) 4*(-1)^n/n^2*cos(n*x)*exp(-n^2*t);
  10.  
  11. % Create a mesh of points between two limits
  12. x0 = pi;
  13. x = linspace(0, x0, 1000);
  14.  
  15. % Create a mesh of points in time
  16. t = linspace(0, 5, 200);
  17.  
  18. figure(1); % Open the figure
  19. plot(x, x.^2, 'b', 'LineWidth', 2); % Plot the base function
  20. ylim([-0.2,pi^2]); % Set the y limits
  21. xlim([0, x0]); % Set the x limits
  22. xlabel('x'); ylabel('u(x,t)');
  23. hold on
  24.  
  25. for j = 1:length(t)
  26. tj = t(j);
  27. u = 1/2*(2*pi^2/3);
  28. for n = 1:N
  29. u = u + R(n, tj, x);
  30. end
  31. if j == 1
  32. p = plot(x, x.^2, 'r');
  33. else
  34. set(p, 'YData', u);
  35. end
  36. drawnow
  37. title(['t = ', num2str(tj)]);
  38. pause(0.1);
  39. if j == 1
  40. pause
  41. end
  42. end