Text and code provided under a Creative Commons Attribution license, CC-BY. (c) Lorena A. Barba, 2013. Thanks: Gilbert Forsyth for help writing the notebooks. NSF for support via CAREER award #1149784. import numpy as np import sympy from sympy import init_printing init_printing(use_latex=True) x, nu, t = sympy.symbols('x nu t') phi = sympy.exp(-(x-4*t)**2/(4*nu*(t+1))) + sympy.exp(-(x-4*t-2*np.pi)**2/(4*nu*(t+1))) phi phiprime = phi.diff(x) phiprime print phiprime from sympy.utilities.lambdify import lambdify u = -2*nu*(phiprime/phi)+4 print u ufunc = lambdify((t, x, nu), u) print ufunc(1,4,3) import matplotlib.pyplot as plt ###variable declarations nx = 101 nt = 100 dx = 2*np.pi/(nx-1) nu = .07 dt = dx*nu x = np.linspace(0, 2*np.pi, nx) #u = np.empty(nx) un = np.empty(nx) t = 0 u = np.asarray([ufunc(t, x0, nu) for x0 in x]) u plt.figure(figsize=(11,7), dpi=100) plt.plot(x,u, marker='o', lw=2) plt.xlim([0,2*np.pi]) plt.ylim([0,10]) for n in range(nt): un[:]=u[:] for i in range(nx-1): u[i] = un[i] - un[i] * dt/dx *(un[i] - un[i-1]) + nu*dt/dx**2*\ (un[i+1]-2*un[i]+un[i-1]) u[-1] = un[-1] - un[-1] * dt/dx * (un[-1] - un[-2]) + nu*dt/dx**2*\ (un[0]-2*un[-1]+un[-2]) u_analytical = np.asarray([ufunc(nt*dt, xi, nu) for xi in x]) plt.figure(figsize=(11,7), dpi=100) plt.plot(x,u, marker='o', lw=2, label='Computational') plt.plot(x, u_analytical, label='Analytical') plt.xlim([0,2*np.pi]) plt.ylim([0,10]) plt.legend() from IPython.core.display import HTML def css_styling(): styles = open("../styles/custom.css", "r").read() return HTML(styles) css_styling()