# -*- coding: utf-8 -*- # 3.0 # # The One Dimensional Heat Equation # from sympy import * from sympy.utilities.lambdify import lambdify init_printing() # # The one dimensional heat conduction # # If $u(x,t)$ is the temperature at position $x$ and time t the one dimensional heat equation is given by: # x = Symbol('x') t = Symbol('t') u = Function('u')(x,t) heat = Eq(u.diff(t),u.diff(x,x)) pprint(heat) # # In the lectures, we have found the particular solution # k = Symbol('k') u_part = exp(-(k*pi)**2*t)*sin(k*pi*x) pprint(u_part) # # if $k$ is an integar, # u_part = u_part.subs(k,1) pprint(u_part) # # we can see that u_part satifies the heat equation # heat = heat.subs(u,u_part) pprint(heat.doit()) # # and the following boundary conditions: # u_0 = u_part.subs(x,0) u_1 = u_part.subs(x,1) print "u_0 = " + pretty(u_0) print "u_1 = " + pretty(u_1) # # and initial conditions: # u_init = u_part.subs(t,0) print "u_init = " + pretty(u_init) # # And here's a 3D-plot of the solution # import numpy as np from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib import cm lam_u_part = lambdify((x,t), u_part) Nx = 10 Nt = 10 xa = np.linspace(0,1,Nx) ta = np.linspace(0,0.5,Nt) X, T = np.meshgrid(xa, ta) Q = [[lam_u_part(Xe,Te) for (Xe, Te) in zip(Xrow,Trow)] for (Xrow,Trow) in zip(X,T)] fig = plt.figure() ax = fig.gca(projection='3d') ax.plot_surface(T, X, Q, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True) fig = plt.figure() # # Numerical solution # # Now, we will use the explicite time stepping scheme derived in the lectures to approximately # compute a solution. # # First, we define the number of unknowns, the mesh size, and the size and number of the time # steps: # n = 20; h = 1./n tau = 1./800; steps = 100 v = np.zeros((n+1,steps+1)) # # Then, we initialize our numerical solution: # lam_u_init = lambdify(x,u_init) for j in range(1,n): v[j,0] = lam_u_init(j*h) # # and set the boundary conditions: # for m in range(0,steps+1): v[0,m] = 0 v[n,m] = 0 # # To compute the numerical solution, we can then use the loop # for j in range(1,n): v[j,1] = v[j,0] + (tau/(h**2)) * ( v[j-1,0] - 2*v[j,0] + v[j+1,0] ) # # We may take a little time to plot the results: # plt.figure() plt.plot(v[:,0]) plt.hold(True) plt.plot(v[:,1]) plt.hold(False) plt.show() # # To compute the solution for all time steps, we use the following loop: # for m in range(0,steps): for j in range(1,n): v[j,m+1] = v[j,m] + tau/h**2 * ( v[j-1,m] - 2*v[j,m] + v[j+1,m] ) # # And plot the results again: # plt.figure() plt.hold(True) for m in range(0,steps+1): plt.plot(v[:,m]) plt.hold(False) plt.show() # # Of course, we would like to compare the numerical solution with the anayltic one. The following # plot shows, that we can come pretty close ... # exact_plot = [lam_u_part(j*h,steps*tau) for j in range(0,n+1)] plt.figure() plt.hold(True) plt.plot(v[:,steps]) plt.plot(exact_plot) plt.hold(False) plt.show()