Context: The atmosphere of Venus wasn’t always as warm as now: a few hundred million years ago it was much more habitable. I am interested in how fast the interior of the planet warmed up, after it’s climate changed.

Model:
Venus is a uniform sphere with some heat conductivity $alpha$, temperature distribution $f(r)$ at $t=0$. I am interested in how the temperature changes ie in $u(r,t)$.

---

Setup

Equation to solve:$$dot{u}=alpha nabla^2 u$$

In spherical polars:
$$dot{u} = alpha frac{1}{r^2} frac{partial}{partial r}left(r^2frac{partial u}{partial r}right)$$

Initial & Boundary conditions:
$$u(r leq R, t=0) = f(r)$$ and $$u(r=R,t geq 0) = f(R)$$
(Ie planet has some initial temperature distribution straight after its outermost layers warmed up, we keep the outermost layers at fixed temperature $f(R)$.)

Change of variables: I change from $u$ to $v(r,t)=u(r,t)-frac{r}{R}f(R)-frac{R-r}{R}f(0)$, so IC satisfied by $v$: $$v(r leq R,t=0)=f(r)-frac{r}{R}f(R) = g(r)$$
BC satisfied by $v$:
$$v(r=R)=0$$
(and $v(r=0)$ will also be $0$ which is ideal).

---

Solution

Strategy: separation of variables. Ie: $v(r,t) = T(t)rho(r)$. Plug this into heat equation, obtain two ODEs:

Temporal ODE

$$T’=-lambda alpha T$$
ie:
$$T = e^{-lambda alpha t}$$
(Will deal with the arbitrary multiplicative constant in the spatial part.)

Spatial ODE

The equation we get after separating variables:

$$rrho”+2rho’+lambda r rho = 0$$

The allowable solution to this is:

$$rho(r)=frac{c_0}{sqrt{lambda}r}sin(sqrt{lambda}r)$$

(The other solution is with $cos$ instad of $sin$, but that goes to infinity as $r$ approaches $0$, so I don’t use those.)

BC requires that:

$$lambda=frac{n^2pi^2}{R^2}$$

IC requires that:

We find a linear combination of $rho_n$s which sum to $g(r)$.

Rearrange spatial equation:

$$-rho”-frac{2}{r}rho’ = lambda rho$$

Recall that Sturm-Liouville operators have orthogonal eigenfunctions, we let’s turn this into SL form by multiplying by $r^2$. Rearrange, get SL form:

$$-r^2rho”-2rrho’=-frac{d}{dr}left(r^2frac{d}{dr}rhoright) = lambda r^2 rho$$

Need normalized eigenfunctions, ie:

$$int_0^R rho(r,n)rho(r,m) r^2 dr = delta_{nm}$$

where $r^2$ has been included because it is the weight function. Arrive at:
$$c_0 = frac{npi}{R}sqrt{frac{2}{R}}$$

To find the linear combination, I use the Green’s function method. $rho(r)$ is given by:

$$rho(r)=int_0^R G(r, zeta)zeta^2rho(zeta) dzeta$$

where $zeta^2$ is the weight function and

$$G(r, zeta) = sum_{n=1}^{infty}frac{1}{lambda_n}rho_n(r)rho_n(zeta)$$

Rearrange to make it clearer that we have a weighted sum, get $rho$:

$$rho(r)=sum_1^{infty}left(frac{1}{lambda_n}int_0^Rrho_n(zeta)zeta^2g(zeta)dzetaright)rho_n(r)$$

Unite results

Bring back $f$, put time dependence back in, transform back to $u$ by adding back what we subtracted from it:

$$u(r,t)=sum_1^{infty}left(frac{1}{lambda_n}int_0^Rrho_n(zeta)zeta^2left(f(zeta)-frac{zeta}{R}f(R)right) d zeta right)
rho_n(r) e^{-lambda_n alpha t} + frac{r}{R}f(R)$$

---

Checking

I would like to check what I have done above is numerically correct, with Python.

# Useful packages
import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm
import scipy.integrate as integrate

Variables:

step=0.001
start=0+step
stop=4
r = np.arange(start, stop+step, step)
r_xsi = np.arange(start, stop+step, step)
R=stop

Defining $lambda_n$, $rho_n$:

def LAMBDA(n,R=stop):
    return n**2*np.pi**2/(R**2)

def rho_func(r,n,R=stop):
    return (n*np.pi/R)*(2/R)**0.5*1/((LAMBDA(n)**0.5)*r)*np.sin((LAMBDA(n)**0.5)*r)

An example initial temperature distribution, $f(r)$:

def f(r,R=stop):
    return (np.exp((20*(r-0.95*R))/R)+np.cos(r/R)**2+1/(r+0.4))

Which transforms into, after changing variables:

def g(r, R=stop, f=f):
    return f(r) - r/R*f(R)

Define the max $n$ we are calculating, calculate $lambda$s:

Nmax=50
ns=range(1, Nmax+1)
LAMBDAS=np.array([LAMBDA(n) for n in ns])

Calculate the integral which will be the weight in the sum:

integrals = np.array([res for res, error in [integrate.quad(lambda zeta: rho_func(zeta, n)*(zeta**2)*g(r=zeta, R=R), 0, R, limit=1000) for n in ns]])

Divide result of $n$th integration by $lambda_n$, multiply by $rho_n$, sum over $n$ just us our formula requires:

res = np.sum(np.stack([coeff*rho_func(r,n) for n, coeff in enumerate(1/LAMBDAS*integrals,1)]),axis=0)

Without taking into account the temporal variation, let’s check if the initial condition (ie when $T(t=0)=1$ anyway) is satisfied:

plt.plot(r,f(r),label='f')
plt.plot(r,res+ r/R*f(R),label='Should be f')
plt.legend()

Which gives:

which is clearly wrong.

---

Question

What am I doing wrong? Are the calculations wrong, or I am just implementing the checking in a wrong way, and the calculations are correct? I am curious.