Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Views: 1334
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
ut=1.142ux2,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)
XX=T1.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 μ=(2n1)π2,n1\mu=\frac{(2n-1)\pi}{2},\quad n\ge 1
Portanto,
Autovalores: λn=((2n1)π2)2\lambda_n=(\frac{(2n-1)\pi}{2})^2
Autofunções: Xn(x)=sin(2n1)π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.14000000000000lnT(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(5750lnt)\displaystyle C e^{\left(-\frac{57}{50} \, l_{n} t\right)}
A solução geral é: u(x,t)=n=1bne1.14lntsin(2n1)π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=2101(x2+2x)sin(2n1)π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π3n312π3n2+6π3nπ3\displaystyle \frac{32}{8 \, \pi^{3} n^{3} - 12 \, \pi^{3} n^{2} + 6 \, \pi^{3} n - \pi^{3}}
π3(2n1)3\displaystyle \pi^{3} {\left(2 \, n - 1\right)}^{3}
32π3(2n1)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︠