Sharednotebooks / ProblemasHomogeneos-SeparacaoVariaveis.sagewsOpen in CoCalc
Download, Raw, Embed
Author: Sergio Luis Lopes Verardi
O Sage como ferramenta auxiliar na aplicação do Método de Separação de Variáveis
Infelizmente "ainda" não é possível simplesmente definir o problema a valores inicial e de fronteira e solicitar ao Sage que o resolva.
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.
E, naturalmente, para construir os gráficos das soluções encontradas.
Resolva o problema homogêneo
∂u∂t=1.14∂2u∂x2,x∈Ω=[0,1],t>0u(0,t)=0,ux(1,t)=0,u(x,0)=−x2+2x\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}
Trata-se de um problema homogêneo que pode ser resolvido por Separação de Variáveis.
u(x,t)=X(x)T(t)u(x,t)=X(x)T(t)
X′′X=T′1.14T=−λ⇒{T′+1.14λT=0X′′+λX=0\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.
É simples verificar que só existem soluções não triviais para λ>0\lambda\gt 0. Então, vamos assumir que λ=μ2\lambda=\mu^2. Logo,
X′′+μ2X=0X''+\mu_2X=0
X(0)=0X(0)=0
X′(1)=0X'(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)
μ2X(x)+D[0,0](X)(x)=0\displaystyle \mu^{2} X\left(x\right) + D[0, 0]\left(X\right)\left(x\right) = 0
K2cos⁡(μx)+K1sin⁡(μx)\displaystyle K_{2} \cos\left(\mu x\right) + K_{1} \sin\left(\mu x\right)
Precisamos, agora, aplicar as condições de contorno.
Aplicando X(0)=0X(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
K2=0\displaystyle K_{2} = 0
(_K1, _K2, mu, x)
K1sin⁡(μx)\displaystyle K_{1} \sin\left(\mu x\right)
Aplicando X′(1)=0X'(1)=0, temos:
eqCC2=diff(solX1,x).substitute(x==1)==0; show(eqCC2)
K1μcos⁡(μ)=0\displaystyle K_{1} \mu \cos\left(\mu\right) = 0
Portanto, uma vez que μ>0\mu > 0, para termos solução não trivial, devemos impor cos⁡μ=0\cos{\mu}=0.
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)
[mu == 1/2*pi]
Ou seja, nós mesmos precisamos resolvê-la; sabemos que μ=(2n−1)π2,n≥1\mu=\frac{(2n-1)\pi}{2},\quad n\ge 1
Portanto,
Autovalores: λn=((2n−1)π2)2\lambda_n=(\frac{(2n-1)\pi}{2})^2
Autofunções: Xn(x)=sin⁡(2n−1)π2xX_n(x)=\sin{\frac{(2n-1)\pi}{2}x}
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())
1.14000000000000 lnT(t)+D[0](T)(t)=0\displaystyle 1.14000000000000 \, l_{n} T\left(t\right) + D[0]\left(T\right)\left(t\right) = 0
Ce(−5750 lnt)\displaystyle C e^{\left(-\frac{57}{50} \, l_{n} t\right)}
A solução geral é: u(x,t)=∑n=1∞bne−1.14lntsin⁡(2n−1)π2xu(x,t)=\displaystyle\sum_{n=1}^{\infty}b_n e^{-1.14 l_n t}\sin{\frac{(2n-1)\pi}{2}x}
O coeficiente bnb_n é calculado impondo-se a condição inicial e usando a propriedade de ortogonalidade das autofunções: bn=21∫01(−x2+2x)sin⁡(2n−1)π2xdxb_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)
328 π3n3−12 π3n2+6 π3n−π3\displaystyle \frac{32}{8 \, \pi^{3} n^{3} - 12 \, \pi^{3} n^{2} + 6 \, \pi^{3} n - \pi^{3}}
π3(2 n−1)3\displaystyle \pi^{3} {\left(2 \, n - 1\right)}^{3}
32π3(2 n−1)3\displaystyle \frac{32}{\pi^{3} {\left(2 \, n - 1\right)}^{3}}
A solução está completa; podemos esboçar seu gráfico.
No primeiro gráfico, a curva em verde é a distribuição inicial de temperatura.
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→∞)=0T(t\rightarrow\infty)=0.
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)
3D rendering not yet implemented
︠b45de380-9645-462d-9184-b2520499e38e︠