27  Problem set 4

Note

Note that Q1 and Q4 are new to complement the coverage on numerical solutions of PDEs.`

Q1. Finite difference formulae

The finite-difference formulae used in our numerical ODE and PDE schemes can be derived using Taylor series approximations. Consider a function \(U\) that can be expanded as a Taylor series.

  1. Show that \[ \begin{gathered} U(x + h) = U(x) + hU'(x) + \frac{h^2}{2}U''(x) + \frac{h^3}{6}U'''(x) + \cdots, \\ U(x - h) = U(x) - hU'(x) + \frac{h^2}{2}U''(x) - \frac{h^3}{6}U'''(x) + \cdots. \end{gathered} \]

  2. By subtracting the two equations above, derive the centred difference formula used for the first-derivatives in our algorithms: \[ U'(x) = \frac{U(x + h) - U(x - h)}{2h} + O(h^2). \] This is what’s known as the centred difference formula for the derivative.

  3. By adding the two Taylor series equations above, derive the centred difference formula used for the second-derivatives in our algorithms: \[ U''(x) = \frac{U(x + h) - 2U(x) + U(x - h)}{h^2} + O(h^2). \]

  4. The centered difference formula for \(U'(x)\) will not immediately work if it is applied to find the derivative at the first mesh point in a discretisation, e.g. \(x = x_0\). Show that the forward difference formula at \(x = x_0\) can be derived as \[ U'(x_0) = \frac{-3U(x_0) + 4U(x_1) - U(x_2)}{2h} + O(h^2). \]

Hint: consider as well the Taylor series for \(U(x + 2h)\) and then combine the expansions for \(U(x + h)\) and \(U(x + 2h)\) to derive the result.

Q2. The wine cellar problem I

\(\nextSection\)

From Chapter 11, we have the heat system: \[ \begin{gathered} \frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2}, \\ T(x = 0, t) = A \cos(\omega t), \\ |T(x, t)| \text{ is bounded as $x \to \infty$}. \end{gathered} \] for \(x \geq 0\). The initial condition can be ignored.

The solution was sought in the form of \(T(x, t) = G(t) H(x)\) using a separation of variables procedure leading to \(\frac{G'}{G} = \kappa \frac{H''}{H} = \lambda\). Go through an argument that the cases of \(lambda = 0\) and \(\lambda\) real lead to incompatible solutions. Review the calculation that \(\lambda = i\omega\) leads to the desired solution given in the lecture notes.

Q3. The wine cellar problem II

\(\nextSection\)

Review the code written on the Noteable repository that implements the explicit Euler scheme for the wine cellar problem.

  • Verify how well the numerical solution matches the exact analytical solution by overlaying the numerical solution with the analytical solution.
  • Investigate the solutions under different choices of parameters and time-stepping tolerances until you are confident you know how the numerical algorithm works.

Q4. von Neumann analysis

\(\nextSection\)

Use von Neumann analysis to determine the stability of the explicit Euler scheme for the 4th order diffusion equation: \[ u_t = -u_{xxxx}. \] You may use the fact that the centred difference approximation of the fourth derivative is: \[ u_{xxxx}(x_j) \approx \frac{u_{j-2} - 4u_{j-1} + 6u_j - 4u_{j+1} + u_{j+2}}{\Delta x^4}. \]

The stability criterion should be \[ \Delta t \leq \frac{\Delta x^4}{8}. \]

Q5. Variable Sun output

\(\nextSection\)

Satellite data indicates that \(Q\), varies roughly between 341.37 W/m2 and 341.75 W/m2, with a period of about 11 years.

  1. Use the simple EBM (Equation 12.3), given by \[ Q(1 - a) = \sigma \gamma T^4, \] with a constant albedo, \(a = 0.3\) and greenhouse gas factor \(\gamma = 0.6\) to estimate the resultant variation (max and min) in the Earth’s mean surface temperature \(T\).

  2. Similar to (a) but this time, use the Budyko balance equation, \[ Q(1 - a) = A + BT \] with \(A = 203.3 \, \mathrm{W} \mathrm{m}^{-2}\) and \(B = 2.09 \, \mathrm{W}/(\mathrm{m}^{2} \, {}^\circ \mathrm{C})\) to estimate the resultant variation in the surface temperature. Use \(a = 0.3\).

  3. The actual variation in surface temperature is in fact less than what you computed above. Why might this be?

Q6. Phase line analysis

\(\nextSection\)

Consider the energy balance equation \[ C \frac{\mathrm{d}T}{\mathrm{d}t} = Q(1 - a(T)) - \sigma \gamma T^4 \equiv G(T). \tag{27.1}\] with \(a\) given by (Equation 12.2). Because the differential equation is autonomous, we can apply phase-line analysis in order to qualitatively understand the evolution. Below is a plot of the function \(G\):

import numpy as np
import matplotlib.pyplot as plt

Q = 342; sigma = 5.67e-8; gam = 0.6;
a = lambda T: 0.5 - 0.2*np.tanh(0.1*(T-265));
T = np.linspace(220, 310, 100);
G = Q*(1-a(T)) - sigma*gam*T**4
fig, ax = plt.subplots()
ax.plot(T, G)
ax.grid(); ax.minorticks_on();
# Customize the major grid
ax.grid(which='major', linestyle='-', linewidth='0.5', color='red')
# Customize the minor grid
ax.grid(which='minor', linestyle=':', linewidth='0.5', color='black')
plt.xlabel('T (in Kelvins)'); plt.ylabel('G(T)');

  1. Sketch the solution \(T(t)\) of this equation for \(t > 0\) if \(T(0) = 230, 240, 260, 270\) and \(300\).

  2. Sketch the solution \(T(t)\) of this equation for \(t > 0\) if \(T(0) = 285\). Then sketch the solution of this equation with the same initial data in the same coordinate system if \(C\) is twice as large. Explain your answer.

  3. If \(\gamma\) is decreased due to the increased greenhouse effect, the entire curve is shifted upwards. Sketch the solution if \(T(0) = 280\). Sketch the solution with the same initial data if \(\gamma\) is decreased. Explain your answer.