ProblemasHomogeneos-SeparacaoVariaveis.ipynb

Kernel: SageMath (stable)

%html
<font size=5>O Sage como ferramenta auxiliar na aplicação do Método de Separação de Variáveis</font><br>
Infelizmente "ainda" não é possível simplesmente definir o problema a valores inicial e de fronteira e solicitar ao Sage que o resolva.<br>
Então, o Sage deve ser utilizado como uma ferramenta auxiliar a ser empregada nas etapas que envolvem cálculos longos ou complicados como a determinação dos coeficientes de Fourier e a solução de EDOs.<br>
E, naturalmente, para construir os gráficos das soluções encontradas.<br>

%html
Resolva o problema homogêneo<br>
$\begin{array}{lcl}
\displaystyle \frac{\partial u}{\partial t}&=&1.14\displaystyle \frac{\partial^2 u}{\partial x^2},\quad x\in\Omega=[0,1],t>0\\
u(0,t)&=&0,\\
u_x(1,t)&=&0,\\
u(x,0)&=& -x^2+2x\\
\end{array}$

%html
Trata-se de um problema homogêneo que pode ser resolvido por Separação de Variáveis.<br>
$u(x,t)=X(x)T(t)$<br>
$\displaystyle\frac{X''}{X}=\displaystyle\frac{T'}{1.14T}=-\lambda\Rightarrow \left\{\begin{array}{l}T'+1.14\lambda T=0\\X''+\lambda X=0\end{array}\right.$<br>
É simples verificar que só existem soluções não triviais para $\lambda\gt 0$. Então, vamos assumir que $\lambda=\mu^2$. Logo,<br>
$X''+\mu_2X=0$<br>
$X(0)=0$<br>
$X'(1)=0$

# Solução da EDO associada a X(x)
mu,x=var('mu','x')
assume(mu>0)
X=function('X')(x)
#mu2=mu^2
edoX=diff(X,x,2)+mu^2*X==0;show(edoX)
solX=desolve(edoX,[X,x]);show(solX)

%html
Precisamos, agora, aplicar as condições de contorno.<br>
Aplicando $X(0)=0$, vem

eqCC1=solX.substitute(x==0)==0; show(eqCC1) # em x=0
solX.variables();
k2=solX.variables()[1]
solX1=solX.substitute(k2==0);show(solX1) #aplica K2=0 na solução

%html
Aplicando $X'(1)=0$, temos:

eqCC2=diff(solX1,x).substitute(x==1)==0; show(eqCC2)

%html
Portanto, uma vez que $\mu > 0$, para termos solução não trivial, devemos impor $\cos{\mu}=0$.<br>
Poderíamos pensar em resolver esta equação, através da função solve() mas, como mostra o exemplo a seguir, o Sage não fornece a solução geral dessa equação trigonométrica:

solve(cos(mu)==0,mu)

%html
Ou seja, nós mesmos precisamos resolvê-la; sabemos que $\mu=\frac{(2n-1)\pi}{2},\quad n\ge 1$<br>
Portanto,<br>
Autovalores: $\lambda_n=(\frac{(2n-1)\pi}{2})^2$<br>
Autofunções: $X_n(x)=\sin{\frac{(2n-1)\pi}{2}x}$

%html
Agora, resolvemos a EDO em T, já levando em conta a existência dos autovalores

# Solução da EDO associada a T(T)
n,t,l_n=var('n','t','l_n')
assume(n,'integer')
assume(n>0)
#l_n=((2*n-1)*pi/2)^2
T=function('T')(t)
edoT=diff(T,t,1)+1.14*l_n*T==0;show(edoT)
solT=desolve(edoT,[T,t]);show(solT.simplify_full())

%html
A solução geral é:
$u(x,t)=\displaystyle\sum_{n=1}^{\infty}b_n e^{-1.14 l_n t}\sin{\frac{(2n-1)\pi}{2}x}$<br>
O coeficiente $b_n$ é calculado impondo-se a condição inicial e usando a propriedade de ortogonalidade das autofunções:
$b_n=\frac{2}{1}\displaystyle\int_0^1 (-x^2+2x)\sin{\frac{(2n-1)\pi}{2}x}dx$

bn1=(2)*integral((-x^2+2*x)*sin((2*n-1)*pi/2*x),v=x,a=0,b=1).simplify_full();show(bn1)
bn2=bn1.denominator().factor();show(bn2)
bn=32/bn2;show(bn)

%html
A solução está completa; podemos esboçar seu gráfico.<br>
No primeiro gráfico, a curva em verde é a distribuição inicial de temperatura.<br>
Como era de se esperar, já que se trata de um fenômeno de difusão e pelas condições de contorno, com o passar do tempo a temperatura tende para o estado estacionário e, neste caso, $T(t\rightarrow\infty)=0$.<br>
Estas propriedades podem ser vistas de uma só vez no gráfico em 3D.

# Calcula a solução u(x,t) com N termos no somatório
xx=var('xx')
N=5
bbn=[32/(((2*n-1)^3)*pi^3) for n in range(1,N+1)]
u=sum(bbn[n]*exp(-1.14*((2*(n+1)-1)/2*pi)^2*t)*sin((2*(n+1)-1)/2*pi*xx) for n in range(0,N))

# Gráfico de u(x,t) para diferentes instantes
P2=plot(u.substitute(t==0),[xx,0,1],color='green');
P3=plot(u.substitute(t==0.01),[xx,0,1]);
P4=plot(u.substitute(t==0.08),[xx,0,1]);
P5=plot(u.substitute(t==0.2),[xx,0,1]);
P6=plot(u.substitute(t==0.6),[xx,0,1]);
P7=plot(u.substitute(t==1.0),[xx,0,1]);
P8=plot(u.substitute(t==1.5),[xx,0,1]);
P9=plot(u.substitute(t==2.0),[xx,0,1]);
(P2+P3+P4+P5+P6+P7+P8+P9).show()

# Gráfico 3D de u(x,t)
plot3d(u,[xx,0,1],[t,0,2],plot_points=[20,20],aspect_ratio=(1,1,0.2),mesh=True)