Sharednotebooks / ProblemasHomogeneos-SeparacaoVariaveis.ipynbOpen in CoCalc
Author: Sergio Luis Lopes Verardi
In [ ]:
%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>
In [ ]:
%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}$
In [ ]:
%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$
In [ ]:
# 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)
In [ ]:
%html Precisamos, agora, aplicar as condições de contorno.<br> Aplicando $X(0)=0$, vem
In [ ]:
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
In [ ]:
%html Aplicando $X'(1)=0$, temos:
In [ ]:
eqCC2=diff(solX1,x).substitute(x==1)==0; show(eqCC2)
In [ ]:
%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:
In [ ]:
solve(cos(mu)==0,mu)
In [ ]:
%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}$
In [ ]:
%html Agora, resolvemos a EDO em T, já levando em conta a existência dos autovalores
In [ ]:
# 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())
In [ ]:
%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$
In [ ]:
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)
In [ ]:
%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.
In [ ]:
# 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)
In [ ]:
In [ ]: