Resolver la ecuación del calor en polares esféricos con condiciones de contorno no homogéneas

Tratando de encontrar la solución en serie de$r\rho''+2\rho'-\lambda r \rho = f$con ciertos IC y BC. Pregunta: qué está mal con la estrategia de solución (la verificación numérica no confirma el resultado).

Mi pregunta:

¿Qué función describe el cambio de temperatura interna de un cuerpo esférico con algún perfil de temperatura inicial cuando su entorno se calienta?

Plan

1. Establezca la ecuación que gobierna la transferencia de calor
2. Ponlo en coordenadas polares esféricas
3. Configurar condiciones iniciales y de contorno
4. Obtenga dos ecuaciones separadas usando la separación de variables, como aquí, en una configuración diferente
5. Obtenga soluciones para aquellos: [AQUÍ ES DONDE VA MAL]
   • Temporal, con suerte, será una EDO de primer orden, algunas exponenciales como solución.
   • La espacial probablemente será más complicada, use: solución en serie
6. Aplique IC y BC para fijar constantes en las soluciones generales, obtenga el resultado final

Ejecución

1. $\dot{u} = \alpha \Delta u$,$R \geq r \geq 0$,$R$siendo el radio de la esfera,$t \geq 0$. Como esperamos una solución esféricamente simétrica,$u=u(r,t)$.$u$no depende de$\theta$&$\phi$.
2. $\dot{u} = \alpha \frac{1}{r^2} \frac{\partial}{\partial r}(r^2\frac{\partial u}{\partial r})$. Desde$u$es sólo una función de$r$&$t$, podemos ignorar la$\theta$&$\phi$dependencia de$\Delta$.
3. • CI:$u(r \leq R,t=0) = f(r)$(es decir, la temperatura interna inicial del objeto esférico)
   • ANTES DE CRISTO:$u(r=R,t>0) = T_f$(es decir, la temperatura del entorno calentado, el calentamiento se toma como instantáneo)
4. Dejar$u=T(t)\rho(r)$. Ecuaciones separadas:
   • temporal:$T' = \lambda \alpha T$
   • espacial:$r\rho''+2\rho'-\lambda r \rho = 0$(similar a esto y esto )
5. Soluciones:
   • Temporal:$T(t)=Ae^{\lambda \alpha t}$
   • espacial: dejar$\rho(r) = \sum_0^{\infty}c_n r^n$. (No permitir poderes negativos porque nos gustaría tener una solución real en$r=0$también.) Sustituya esto en la ecuación para$\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$$ que es igual a: $$\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$$ En la primera suma, cambia cada$n$a$n+2$& baje el límite inferior de la suma en 2. En la segunda suma, cambie cada$n$a$n+1$& baje el límite inferior de la suma en 1. Obtenga: $$\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$$ La ecuación anterior solo puede ser verdadera si es verdadera para cada potencia de$r$. Considera el$r^{n+1}$términos: $$(n+2)(n+1)c_{n+2} - \lambda c_n = -2(n+2)c_{n+2}$$ Reorganizar: $$(n+4)(n+1)c_{n+2}=\lambda c_n$$ Volver a escribir: $$(n+2)(n-1)c_n = \lambda c_{n-2}$$ Obtener la relación de recurrencia: $$c_n = \frac{\lambda c_{n-2}}{(n+2)(n-1)}$$ Para no meterse en problemas con la división cero, establezcamos$c_{-1}$a$0$, renderizando todo$c_{n=odd}$cero. Somos libres de elegir$c_0$, determinando todos los coeficientes pares. ¡Tracemos esto y comprobemos que estamos en lo cierto! Usando$c_0=\lambda=1$, calculó una aproximación a$\rho$usando$c_n$depende de$n=100$. También calculó una aproximación a$r\rho''+2\rho'-\lambda r \rho$, que debe ser cero. (Cuaderno disponible aquí , el código también se encuentra en la parte inferior de esta publicación). Desafortunadamente, no lo es:

(Tenga en cuenta que esta publicación tiene una ecuación homogénea muy similar a la mía, pero no usa soluciones en serie, pero afirma que esta ecuación "se integra en" la solución).

6. Obtenemos:$u=T(t,A)\rho(c_0,r)$. Incluso si mi$\rho$la función era correcta, BC apenas parece satisfactoria para el general$f$con estas constantes. ¿Me estoy perdiendo de algo? Recuerdo que en la escuela (/uni) resolvimos el caso homogéneo, cuando$f$es$0$, obtuvo una solución general, luego cambió$f$a algo distinto de cero, encontró una solución particular, luego agregó las soluciones general y particular. Mi plan para esto:

   • Encuentra funciones propias ortonormales$y_n$& valores propios$\sigma_n$para$Ly_n=\sigma_n y_n$, con$L = \frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r}$

   • Arreglar$T(t=0)=1$, reescribe IC como$rL\rho=r(L\rho_g + L\rho_p) = rL\rho_p = f$, es decir$L\rho_p=f/r$dónde$\rho_g$&$\rho_p$son las soluciones generales y particulares en el$u$, dónde$u=T(t)\rho(r)=T(t)(\rho_{general}(r)+\rho_{particular}(r))$,$u$resolviendo$\dot{u} = \alpha \Delta u$.

   • Dejar$\rho_p=\sum_0^{\infty}s_n y_n$

   • Desde

     $$\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}},$$ dónde$\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$. Por lo tanto, puedo calcular todo$s_n$s, por lo tanto obtener$\rho_p$, dándome un$\rho=\rho_g+\rho_p$, con el que puedo satisfacer el BC & IC:$\rho_g$satisfará BC, mientras que$\rho_p$satisfará IC.

Pregunta

¿Qué sale mal en el quinto paso del plan? ¿Es correcto el esquema del último paso, el sexto paso, o me estoy perdiendo algo?

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()