Planetary warmth

Some fannying about with the heat equation

Lets say we have a sphere with some radius-dependent initial temperature distribution inside it. We keep the temperature on the surface of this sphere constant. We would like to know how temperature changes inside. By solving this equation we get closer to answering some questions, such as: how fast can we heat up an asteroid if we don't want to drill much, how a frozen bit of ice melts, and how long can a spherical polar bear last underwater? Questions we are all worried about. (Btw, the main point is that we get better in solving DEs.)

The equation governing temperature distribution

$$\dot{u}=\alpha \nabla^2 u$$ Lets put this into spherical polar coordinates: $$\dot{u} = \alpha \frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2\frac{\partial u}{\partial r}\right)$$ (We care only about $r$ dependence.)

Boundary and Initial Conditions

The initial temperature distribution is our initial condition: $$u(r \leq R, t=0) = f(r)$$ The boundary condition is the temperature outside the body (ie sea around polar bear): $$u(r=R,t \geq 0) = f(R)$$ I am not fond of non-zero boundary conditions. Lets change $u$ to $v$ in a way which make the BCs for $v$ zero: $$v(r,t)=u(r,t)-f(r=R)$$ Then, $v$ will satisfy: $$v(R, t \geq 0) = u(R, t \geq 0)-f(R) = 0$$ We solve for $v$ and leater transform back to $u$. $$v(r \leq R,t = 0) = f(r)-f(r=R) = g(r)$$ $$v(r=R, t \geq 0) = 0$$

Sep Var

Let $u(r,t) = \rho(r)T(t)$. Separated equations we arrive at:

Temporal: $$T'=-\lambda \alpha T$$ Spatial: $$r\rho''+2\rho'+\lambda r \rho = 0$$ $\lambda$ is a constant we don't know yet but we count on BCs to get it.

Temporal ODE

Just observe $T$ to be: $T(t)=Ae^{-\lambda \alpha t}$ for arbitrary $A$. We'll deal with arbirary constant mulitplicative factors in the spatial solution, so let's just drop the $A$: $$T(t)=e^{-\lambda \alpha t}$$

Spatial ODE

We are solving $r\rho''+2\rho'+\lambda r \rho = 0$. Lets do series solution, let $\rho(r) = \sum_0^{\infty}c_nr^n$. We don't want negative powers of $r$ because we want to have a finite $\rho$ at $r=0$. Subsitute in, compare coefficients, stay physical. We arrive at: $$\rho(r)=\frac{c_0}{\sqrt{\lambda}r}\sin(\sqrt{\lambda}r)$$ Our new wonderful BC tells us that $v(R,t \geq 0) = 0$. Since the temporal part is not going to do it (it is an exponential, not even dependent on $r$), the spatial part has to make it $0$. When is $\rho(r=R)$ $0$? When the $\sin$ part is $0$. This is the case when $\sqrt{\lambda}R=n\pi$, ie when $$\lambda=\frac{n^2\pi^2}{R^2}$$

Orthonormality is our friend

Rearrange spatial ODE to: $$-\rho''-\frac{2}{r}\rho' = \lambda \rho$$ LHS is not in Sturm-Liouville form. The weight function which can make it SL type is $w=r^2$. The eigenfunctions of the SL operator are still: $$\rho_n(r)=\frac{c_0}{\sqrt{\lambda_n}r}\sin(\sqrt{\lambda_n}r)$$ These are now orthogonal. To make them orthonormal, we need to find $c_0$. What we want is: $$\int_0^R \rho_n(r)\rho_m(r) w(r) dr = \delta_{nm}$$ Are we integrating over all space or over the all possible values of $r$? The latter, that's why we don't include an additional $r^2\sin(\theta)$ term in the integrand. Think about it deeply, why we don't integrate over all space within planet/droplet/polarbear? We arrive at: $$c_0 = \frac{n\pi}{R}\sqrt{\frac{2}{R}}$$

How much we want from each eigenfunction

Note that all of the eigenfunctions are solutions to the spatial ODE. We need to figure out how much we want from each eigenfunction. All the eigenfunctions are nonzero at $r=0$, and they are all $0$ at $R$; there are no other eigenfunctions (we don't use $\cos$ instead of $\sin$ in them, even though they are solutions too, because we want to stay finite at $r=0$). So the $\rho_n$ eigenfunctions better span the whole function-space. With a bit of notation abuse, we can write: $$g=\sum_{n=1}^{\infty}\langle g | \rho_n \rangle_w | \rho_n \rangle$$ So the amount of we need from each $\rho_n$ is $\langle g | \rho_n \rangle_w$, ie $\int_0^R g(\zeta) \rho_n(\zeta) \zeta^2 d \zeta$.


Combine spatial result with temporal dependence, transform back from $v$ to $u$: $$u(r,t)=\sum_0^{\infty}\int_0^R g(\zeta) \rho_n(\zeta) \zeta^2 d \zeta \rho_n(r) e^{-\frac{n^2\pi^2}{R^2}\alpha t}+f(R)$$


Absolutely not numpy optimized code can be found here. Temperature evolution with an arbitrarily chosen initial temperature distribution, produced by the above code: