34  Problem set 4 solutions

Note

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

Q1. Finite difference formulae

Solution a)

1. Show that $$\begin{array}{c}U(x+h)=U(x)+h{U}^{\prime }(x)+\frac{h^2}{2}{U}^{″}(x)+\frac{h^3}{6}{U}^{‴}(x)+\cdots ,\\ U(x-h)=U(x)-h{U}^{\prime }(x)+\frac{h^2}{2}{U}^{″}(x)-\frac{h^3}{6}{U}^{‴}(x)+\cdots .\end{array}$$

This is straightforwards and by definition of Taylor’s theorem.

Solution b)

2. By subtracting the two equations above, derive the centred difference formula used for the first-derivatives in our algorithms: $$U^{\prime }(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.

Again this is immediate by subtracting the given formulae.

Solution c)

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).$$

This part is also immediate. Adding the two equations gives $$U(x+h)+U(x-h)=2U(x)+h^2{U}^{″}(x)+O(h^3),$$ and simplify.

Solution d)

4. The centered difference formula for $U^{\prime }(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^{\prime }(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.

This is a bit more involved. Let us explain firstly why the standard approach doesn’t work: $$\begin{array}{c}U(x+2h)=U(x)+2h{U}^{\prime }(x)+2h^2{U}^{″}(x)+\frac{4}{3}{h}^3{U}^{‴}(x)+\cdots \\ U(x+h)=U(x)+h{U}^{\prime }(x)+\frac{h^2}{2}{U}^{″}(x)+\frac{h^3}{6}{U}^{‴}(x)+\cdots .\end{array}$$

Now adding gives $$U(x+2h)+U(x+h)=2U(x)+3h{U}^{\prime }(x)+\frac{5h^2}{2}{U}^{″}(x)+\dots$$

However, if we isolate $U^{\prime }(x)$ we obtain a formula where the error term is $O(h{U}^{″}(x))$ instead of $O(h^2)$.

The solution is to consider a linear combination: $$U^{\prime }(x)=aU(x)+bU(x+h)+cU(x+2h),$$ and insert the above formulae. Keeping track of the term multiplying the second derivative, $U^{″}(x)$, we then need: $$b/2+2c=0⟹b=-4c.$$ The remaining two equations are: $$a+b+c=0,bh+2ch=1.$$

The second yields $c=-1/(2h)$ and hence $b=2/h$. The first equation then gives $a=-3/(2h)$.

Q2. The wine cellar problem I

$$

Solution

Separation of variables leads to the equation $$\frac{G^{\prime }}{G}=\kappa \frac{H^{″}}{H}=\lambda .$$

Crucially in considering possible solutions, we want solutions to obey the necessary boundary condition of $T(0,t)=A\cos (\omega t)=G(t)H(0)$. This is highly restrictive.

We first consider $\lambda =0$. Then we have that $G(t)=\text{const.}$. This is then no way to satisfy the time-varying boundary condition.

We next consider $\lambda \in \mathbb{R}$. Then $G(t)=C{e}^{\lambda t}$. So again the behaviour of the time-dependent part is incompatible with the oscillating boundary conditions.

Therefore the only remaining possibility is $\lambda =i\alpha$ where $\alpha$ is real. Then we have that $$G(t)=C_1{e}^{i\alpha t}+C_2{e}^{-i\alpha t}.$$ Written in terms of real-valued functions, this gives essentially a linear combination of $\cos (\alpha t)$ and $\sin (\alpha t)$. In order to be compatible with the boundary conditions, we would thus need $$G(t)=\text{const.}\cos (\omega t)=\text{const.}\times \mathrm{Re}(e^{i\omega t}),$$ i.e. taking the imaginary part (if we wish to continue manipulating exponentials).

Now solving for $H(x)$, we set $H(x)=e^{rx}$ as a solution ansatz, giving the characteristic equation $$r^2=i\omega ⟹r=\pm \sqrt{i\omega }.$$ Therefore the general solution for $H$ yields $$H(x)=D_1{e}^{\sqrt{i\omega }x}+D_2{e}^{-\sqrt{i\omega }x}.$$ In order to parse the above form, it is best to re-write the real and complex parts of the exponential argument separately. We note that by Euler’s theorem, $$\sqrt{i}=(e^{\pi i/2}{)}^{1/2}=e^{\pi i/4}=\frac{1+i}{\sqrt{2}}.$$ So now we can re-write $$H(x)=D_1{e}^{\sqrt{\omega /2}x}{e}^{i\sqrt{\omega /2\kappa }x}+D_2{e}^{-\sqrt{\omega /2\kappa }x}{e}^{-i\sqrt{\omega /2}x}.$$ By the boundedness assumption, we require that $H(x)$ does not blow up as $x\to \mathrm{\infty }$. So therefore $D_1=0$. We can finally put the results together, yielding: $$T(x,t)=\text{const.}\times \mathrm{Re}\left[e^{i\omega t}{e}^{-\sqrt{\omega /2\kappa }x}{e}^{-i\sqrt{\omega /2\kappa }x}\right]=\text{const.}\times \mathrm{Re}\left[e^{-\sqrt{\omega /2\kappa }x}{e}^{i(-\sqrt{\omega /2\kappa }x+\omega t)}\right].$$ Note that the reason you have re-written it with the exponentials arranged in this fashion is because the portion of the exponential with the imaginary argument becomes oscillatory. Taking the real part yields $$T(x,t)=A{e}^{-\sqrt{\omega /2}x}\cos (-\sqrt{\omega /2\kappa }x+\omega t).$$ You can double check that the above solution satisfies all the required conditions, including the condition at $x=0$ (forcing us to choose the prefactor $A$).

Q3. The wine cellar problem II

$$

Solution

Solution script can be found on the solutions/ folder of the Noteable Python website.

The numerical solution is tricky for two reasons. First, you must consider the right range of parameters. Second, the numerical algorithm suffers from instability, so you must consider the right spatial and time discretisation.

To begin, from the notes Chapter 11, you are given that typical parameters are $$\omega =\frac{2\pi }{3.15\times {10}^7\mathrm{s}}.$$ This has already been converted to seconds, so instead use $$\omega =2\pi {\mathrm{yrs}}^{-1}.$$ Indeed this value of $\omega$ was chosen so that $\cos (\omega t)$ has a wavelength of 1 years.

Next, it is given that $$\kappa =0.002{\mathrm{cm}}^2/\mathrm{s}.$$ We use the fact that $$1\mathrm{yr}=3.15\times {10}^7\mathrm{s},1\mathrm{cm}={10}^{-2}\mathrm{m}.$$ So $$\kappa =\frac{2\cdot {10}^{-3}\cdot {10}^{-4}{\mathrm{m}}^2}{3.15\times {10}^{-7}\mathrm{yr}}=(2/3.15){\mathrm{m}}^2{\mathrm{yr}}^{-1}\approx 0.635{\mathrm{m}}^2{\mathrm{yr}}^{-1}.$$ When distributing the typical mesh in the $x$ direction, we want to make sure we go deep enough. We had found that putting the cellar in about 4m deep was enough, so as long as your code solve for a range of $x$ values, say, between $x=0$ and $x=10$, this should be sufficient.

The script is given in the solutions folder.

Q4. von Neumann analysis

$$

This note should be filled in soon. The analysis is similar to the analysis for the heat equation as done in the problem class.

Q5. EBMs and variable sun output

$$

1. Regarding the resultant variation on the Earth’s mean surface.

Solution

The variation seems to be very small. Ignoring units, we have $$T={\left(\frac{1-a}{\sigma \gamma }\right)}^{1/4}{Q}^{1/4}\approx 67.3505\times Q^{1/4}.$$ Substituting the values of $Q$ given this gives a temperature that ranges from $T=289.499\mathrm{K}$ to $289.580\mathrm{K}$. Converting this to Celsius gives $16.35{}^{\circ }\mathrm{C}$ to $16.43{}^{\circ }\mathrm{C}$ i.e. a fraction of a degree.

2. Regarding the variation using the Budyko model.

Solution

We obtain the range of approximately $T=17.06{}^{\circ }\mathrm{C}$ to $17.19{}^{\circ }\mathrm{C}$.

3. Regarding why the disparity between actual surface measurements.

Solution

There are many ways to respond, but one major factor is the oceans, which are ignored in our model. The oceans provide a massive energy sink for the planet, and so we would expect that the minor variations of $Q$ due to the 11-year cycle would be smoothed out by the effects of the oceans.

Q6. Phase line analysis

$$

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

Solution

You should be able to do this question by hand, but the following graph is generated via the accompanying Jupyter script in the solutions folder. The point is that once the steady-state solutions are known (dashed) then each solution given by the different initial condition can be approximated by whether it tends towards or away from the nearest fixed point.

Figure 34.1: Graphical solution for Q3a

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.

Solution

The below numerical solution shows blue for $C=1$ and orange for $C=4$. Notice that increasing $C$ seems to decrease the rate of change of the evolution. Indeed, since $$\frac{\partial T}{\partial t}\propto \frac{1}{C},$$ then multiplying $C$ by factor is equivalent to slowing down the evolution by the inverse of that factor. A factor of $C$ that is twice as large would slow down the evolution by half. Note that the oscillations shown in the graph below are numerical artifacts due to the tolerances on the ODE solver (how do you know this?)

Graphical solution for Q3b

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.

Solution

Decreasing $\gamma$ has the effect of shifting the curve upwards. The hottest steady state $T_3$ is consequently increased. There are two main effects. First, the steady-state is much hotter, so the system will tend towards a much hotter state. Second, because decreasing $ $ will increase the rate of change of $T$ (since we are subtracting less via the factor $-\gamma T^4$), then the evolution towards the hotter state is initially much faster.

You see that in the following diagram.

Graphical solution for Q3c

The above makes sense on a physical level. Decreasing $\gamma$ is equivalent to increasing greenhouse gases. We do expect the system to tend to a hotter state, then.