1834.sagews

odd=lambda f: lambda x: f(x) if x>=0 else -f(-x)
extension=lambda f,t: lambda x: f(float(x+t)%float(2*t)-t)
odd_extension=lambda f,t: extension(odd(f),t)
dalambert=lambda c,L,f,g: lambda x,t: (odd_extension(f,L)(x-c*t)+odd_extension(f,L)(x+c*t))/2+1/(2*c)*numerical_integral(odd_extension(g,L),x-c*t,x+c*t)[0]

def sum2(f,iter):
    ret=0
    for i in iter:
        ret+=f(i).n()
    return ret
a=lambda f,L: lambda n: 2/L*numerical_integral(lambda x: f(x)*cos(n*pi/L*x),0,L)[0]
b=lambda g,L,c: lambda n: 2/(c*n*pi)*numerical_integral(lambda x: g(x)*sin(n*pi/L*x),0,L)[0]
wave_sep=lambda c,L,f,g,N: lambda x,t: sum2(lambda n: sin(n*pi/L*x)*(a(f,L)(n)*cos(c*n*pi/L*t)+b(g,L,c)(n)*sin(c*n*pi/L*t)),range(1,N+1))

f=lambda x: 4*x if x in RealInterval(0,.25) else 4*(1/2-x) if (x in RealInterval(.25,.5)) else 0
g(x)=0

plot3d(wave_sep(1,1,f,g,10),[0,1],[0,1],adaptive=True,mesh=True,dots=True)

plot3d(dalambert(1,1,f,g),[0,1],[0,1],adaptive=True,mesh=True,dots=True)

for t in range(11):
    plot([lambda x: dalambert(1,1,f,g)(x,t*.1),lambda x: wave_sep(1,1,f,g,5)(x,t*.1)],[0,1])

def fd(c,L,f,h,k,max_t):
    s=c^2*k/h^2
    n=floor(L/h)
    max_j=floor(max_t/k)
    mat=[]
    for j in range(max_j):
        mat.append([0]*(n+1))
    for i in range(n+1):
        mat[0][i]=f(i*h)
    for j in range(1,max_j):
        for i in range(1,n):
            mat[j][i]=(1-2*s)*mat[j-1][i]+s*(mat[j-1][i+1]+mat[j-1][i-1])
    return(mat)

b2=lambda f,L: lambda n: 2/L*numerical_integral(lambda x: f(x)*sin(n*pi/L*x),0,L)[0]
heat_sep=lambda c,L,f,N: lambda x,t: sum2(lambda n: exp(-(c*n*pi/L)^2*t)*b2(f,L)(n)*sin(n*pi/L*x),range(1,N+1))

f2=lambda x: 100*sin(3*pi*x)^2 if x in RealInterval(0,1/3) else 45*(1/3-x)

plot3d(heat_sep(1/2,1,f2,5),[0,1],[0,1/4])

list_plot3d(fd(1/2,1,f2,1/20,1/200,1/4),interpolation_type='nn')

plot(f2,0,1)

for curve in fd(1/2,1,f2,1/20,1/200,1/4):
    list_plot(curve)