Grafico-ProblemaNaoHomogeneo1.sagews

t,x=var('t','x')
gam,A,B=var('gam','A','B')
k,L=var('k','L')
gam=15;A=1;B=2;k=1;L=1;
# Para simplificar, condição inicial nula f(x)=0
v=function('v')(x)
v=-(1/2)*(gam/k)*x^2+(1/L)*(B-A+gam*L^2/(2*k))*x+A
# Gráfico da função v(x)
P1=plot(v,[x,0,L],color='red',linestyle='--')
# Número de Termos na Soma Parcial
N=5

# Calcula N coeficientes de Fourier (neste caso, série apenas em senos)
# Calcula usando funções do Sage
#fc = Piecewise([[(0,L),-v]])
#bn=[fc.sine_series_coefficient(n,L) for n in range(1,N+1)]
# Ou calcula usando a fórmula que obtivemos no tutorial teórico
bbn=[(2*gam)/(n^3*pi^3)*((-1)^n-1)-2/(n*pi)*(A-B*(-1)^n) for n in range(1,N+1)]

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

# Gráfico de u(x,t) para diferentes instantes
P2=plot(u.substitute(t==0),[x,0,L]);
P3=plot(u.substitute(t==0.01),[x,0,L]);
P4=plot(u.substitute(t==0.08),[x,0,L]);
P5=plot(u.substitute(t==0.2),[x,0,L]);
P6=plot(u.substitute(t==0.6),[x,0,L],color='green');
(P1+P2+P3+P4+P5+P6).show()

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

#Animação bem simples
P=[]
P.append(plot(u.substitute(t==0),[x,0,L]))
for i in range(1,20):
    intt=0.03*i
    P.append(plot(u.substitute(t==intt),[x,0,L]))
#uanim=P
a=animate(P)
a.show(delay=50)

/ext/sage/sage-8.0/local/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment. warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

3D rendering not yet implemented

t,x=var('t','x')
k,Em,c=var('k','Em','c')
k=1;c=1;Em=100;
v=function('v')(x,t)
v=Em*cos(k*(x-c*t))
P1=plot(v.substitute(t==0),[x,0,pi],color='red',linestyle='--')
P2=plot(v.substitute(t==1),[x,0,2*pi],color='green');
P3=plot(v.substitute(t==2),[x,0,3*pi],color='blue');
P4=plot(v.substitute(t==3),[x,0,4*pi],color='cyan');
#P3=plot(u.substitute(t==0.01),[x,0,L]);
#P4=plot(u.substitute(t==0.08),[x,0,L]);
#P5=plot(u.substitute(t==0.2),[x,0,L]);
#P6=plot(u.substitute(t==0.6),[x,0,L],color='green');
(P1+P2+P3+P4).show()
#(P1+P2+P3+P4+P5+P6).show()

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

#Animação bem simples
#P=[]
#P.append(plot(u.substitute(t==0),[x,0,4*pi]))
#for i in range(1,20):
#    intt=0.03*i
#    P.append(plot(v.substitute(t==intt),[x,0,4*pi]))
##uanim=P
#a=animate(P)
#a.show(delay=50)