Trying to find the series solution of $r\rho''+2\rho'-\lambda r \rho = f$ with certain ICs & BCs. Question: what is wrong with solution strategy (numerical check does not confirm result).
My question:
What function describes the internal temperature change of a spherical body with some initial temperature profile when its surroundings is heated up?
Plan
- Set up equation governing heat transfer
- Put it into spherical polar coordinates
- Set up Initial & Boundary Conditions
- Obtain two separated equations using separation of variables, as
here, in a different setup - Obtain solutions for those: [THIS IS WHERE IT GOES WRONG]
- Temporal one hopefully be a first order ODE, some exponential as solution
- Spatial one probably will be more complicated, use: series solution
- Apply IC & BC to fix constants in the general solutions, get final result
Execution
- $\dot{u} = \alpha \Delta u$, $R \geq r \geq 0$, $R$ being the radius of
the sphere, $t \geq 0$. Since we expect a spherically symmetric
solution, $u=u(r,t)$. $u$ it does not depend on $\theta$ & $\phi$. - $\dot{u} = \alpha \frac{1}{r^2} \frac{\partial}{\partial r}(r^2\frac{\partial u}{\partial r})$. Since
$u$ is only a function of $r$ & $t$, we can ignore the $\theta$ & $\phi$ dependence of $\Delta$. -
- IC: $u(r \leq R,t=0) = f(r)$ (ie initial internal temperature of the spherical object)
- BC: $u(r=R,t>0) = T_f$ (ie temperature of heated up surroundings
, heating up taken as instantaneous)
- Let $u=T(t)\rho(r)$.
Separated equations: - Solutions:
- Temporal: $T(t)=Ae^{\lambda \alpha t}$
- Spatial:
let $\rho(r) = \sum_0^{\infty}c_n r^n$. (Not allowing negative powers because we would like to have a real solution at $r=0$ as well.) Substitute this in to the equation for $\rho$:
$$r\sum_2^{\infty}n(n-1)c_nr^{n-2}+2\sum_1^{\infty}nc_nr^{n-1}-\lambda r\sum_0^{\infty}c_n r^n=0$$
which is equal to:
$$\sum_2^{\infty}n(n-1)c_nr^{n-1}+2\sum_1^{\infty}nc_nr^{n-1}-\lambda \sum_0^{\infty}c_n r^{n+1}=0$$
In the first sum, change every $n$ to $n+2$ & lower the lower limit of summation by 2. In the second sum, change every $n$ to $n+1$ & lower the lower limit of summation by 1. Obtain:
$$\sum_0^{\infty}(n+2)(n+1)c_{n+2}r^{n+1}+2(n+1)c_{n+1}r^{n}-\lambda c_{n}r^{n+1} = 0$$
The above equation can only be true if it is true for each powers of $r$. Consider the $r^{n+1}$ terms: $$(n+2)(n+1)c_{n+2} – \lambda c_n = -2(n+2)c_{n+2}$$
Rearrange:
$$(n+4)(n+1)c_{n+2}=\lambda c_n$$
Rewrite:
$$(n+2)(n-1)c_n = \lambda c_{n-2}$$
Obtain recurrence relation:
$$c_n = \frac{\lambda c_{n-2}}{(n+2)(n-1)}$$
To not to get into trouble with zero division, lets set $c_{-1}$ to $0$, rendering all $c_{n=odd}$ zero. We are free to choose $c_0$, determining all even coefficients.
Let's plot this & check that we are right!
Using $c_0=\lambda=1$, computed an approximation to $\rho$ using $c_n$s up to $n=100$. Also computed an approximation to $r\rho''+2\rho'-\lambda r \rho$, which should be zero. (Notebook available here, the code is also on the bottom of this post.)
Unfortunately, it isn't:
(Note that this post has a very similar homogeneous equation to mine, but it doesn't use series solutions, but claims that this equation "integrates into" the solution.)
-
We obtain: $u=T(t,A)\rho(c_0,r)$. Even if my $\rho$ function was
correct, BC hardly seem satisfiably for general $f$ with these
constants. Am I missing something? I recall in school (/uni) we
solved the homogeneous case, when $f$ is $0$, obtained a general
solution, then changed $f$ to something nonzero, found a particular
solution, then added the general & the particular solutions. My plan for this:-
Find orthonormal eigenfunctions $y_n$ & eigenvalues $\sigma_n$ for $Ly_n=\sigma_n y_n$, with $L = \frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r}$
-
Fix $T(t=0)=1$, rewrite IC as $rL\rho=r(L\rho_g + L\rho_p) = rL\rho_p = f$, ie $L\rho_p=f/r$ where $\rho_g$ & $\rho_p$ are the general & particular solutions in the $u$, where $u=T(t)\rho(r)=T(t)(\rho_{general}(r)+\rho_{particular}(r))$, $u$ solving $\dot{u} = \alpha \Delta u$.
-
Let $\rho_p=\sum_0^{\infty}s_n y_n$
-
Since $$\langle L\rho_p | y_m \rangle = \langle L\sum_{n=0}^{\infty}s_n y_n | y_m \rangle = \langle \sum_{n=0}^{\infty}s_n L y_n | y_m \rangle = \langle \sum_{n=0}^{\infty}s_n \sigma_n y_n | y_m \rangle = s_n \sigma_n \langle \sum_{n=0}^{\infty}y_n | y_m \rangle = s_m \sigma_m
\implies {s_m=\frac{\langle L\rho_p | y_m \rangle}{\sigma_m}},$$ where $\langle f | g \rangle = \int_{all space} f g dV = \int_r=0^{\infty} \int_{\theta=0}^{2\pi} \int_{\phi=0}^{\pi} f(r, \phi, \theta) g(r, \theta, \phi) r^2 \sin(\phi) dr d\phi d\theta$. Therefore, I can compute all $s_n$s, therefore obtain $\rho_p$, giving me a $\rho=\rho_g+\rho_p$, with which I can satisfy the BC & IC: $\rho_g$ will satisfy BC, while $\rho_p$ will satisfy IC.
-
Question
What goes wrong in the 5th step of the plan and is the outline for the last, 6th step is correct, or am I missing something?
The code (Notebook here)
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
LAMBDA = 1
c1=1
def Cn(n, LAMBDA = LAMBDA, c1=c1):
if n%2==1:
return 0
elif n==0:
return 1
else:
return LAMBDA*Cn(n-2)/ ((n+2)*(n-1))
start=-1
stop=10
step=0.0001
rs = np.arange(start, stop+step, step)
rho = [sum([Cn(n)*(r**n) for n in range(0,100+1)]) for r in tqdm(rs)]
def deriv(arr, dx=step):
return np.gradient(arr, dx)
rhoderiv = deriv(rho)
rhoderivderiv = deriv(rhoderiv)
this_is_rather_be_zero = [r*rhoderivderiv[i] + 2*rhoderiv[i] - LAMBDA * r * rho[i] for i, r in enumerate(rs)]
plt.plot(rs, rho,c='k',label='rho')
plt.plot(rs, this_is_rather_be_zero,c='g',label='shouldBeZero')
plt.axvline(c='r')
plt.axhline(c='r')
plt.xlim([-1,5])
plt.ylim([-1,5])
plt.legend()
Best Answer
Firstly I should apologise for being so careless in the comments. I sent us on a silly hunt for coefficients when you were right the first time. But your code did not implement it correctly (see bottom of answer).
In this answer I only deal with the solving of the ODE in step 5, and the issue of the plot showing that it doesn't solve the equation. It may be a good idea to streamline step 6 into a separate question.
From $$ r\rho''+2\rho'-\lambda r \rho = 0$$ Substituting in $\rho=\sum_0^\infty c_n r^n$ we get $$ r\sum_2^\infty n(n-1)c_n r^{n-2}+2\sum_1^\infty n c_nr^{n-1}-\lambda r \sum_0^\infty c_n r^n = 0$$ i.e. $$ \sum_2^\infty n(n-1)c_n r^{n-1}+2\sum_1^\infty n c_nr^{n-1}-\lambda \sum_0^\infty c_n r^{n+1} = 0$$ bumping indices by $\pm1$, and defining $c_{-1} := 0$, $$\sum_1^\infty (n+1)nc_{n+1} r^{n}+2\sum_0^\infty (n+1) c_{n+1}r^{n}-\lambda \sum_0^\infty c_{n-1} r^{n} = 0$$ Finally, we write $\sum_1^\infty (n+1)nc_{n+1} r^{n} = \sum_0^\infty n(n+1)c_{n+1} r^{n}$ and compare coefficients: $$(n+1)n c_{n+1} + 2(n+1) c_{n+1} - \lambda c_{n-1} = 0 $$ i.e. $$(n+1)(n+2) c_{n+1} - \lambda c_{n-1} = 0. $$ Given initial data $c_{-1}=0$, $c_0$ arbitrary we get $c_{\text{odd}}=0$, and setting $n=2k+1$, $$ c_{2(k+1)} = \frac{\lambda c_{2k}}{(2k+2)(2k+3)}$$ Staring hard/Wolfram|Alpha imforms me that the solution is $$ c_{2k} = \frac{\lambda^{k}c_0}{\Gamma(2k+2)}= \frac{\lambda^{k}c_0}{(2k+1)!}, $$ so if $\lambda>0$, $$\rho(r) = \sum_{k=0}^\infty \frac{\lambda^{k}c_0}{(2k+1)!}r^{2k}=c_0\sum_{k=0}^\infty\frac{(\sqrt{\lambda} r)^{2k}}{(2k+1)!}=\frac{c_0}{\sqrt {\lambda} r} \sinh(\sqrt\lambda r)$$
And this is indeed a solution. If $\lambda<0$ then a solution is $\frac1{\sqrt {-\lambda} r} \sin(\sqrt{-\lambda} r)$. If you allow solutions with $c_{-1}\neq 0$, you get also e.g. $\frac{c_0}{\sqrt {\lambda} r} \cosh(\sqrt\lambda r)$.
So what went wrong in the plot? You put in the formula for $c_n$ wrong. It should be
instead of
return LAMBDA*Cn(n-2)/ ((n+2)*(n-1)). This gives e.g.Cn(2)as 1/6 as it should, and the plot