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. from mpl_toolkits.mplot3d import Axes3D ##New Library required for projected 3d plots import numpy as np import matplotlib.pyplot as plt ###variable declarations nx = 81 ny = 81 nt = 100 c = 1 dx = 2.0/(nx-1) dy = 2.0/(ny-1) sigma = .2 dt = sigma*dx x = np.linspace(0,2,nx) y = np.linspace(0,2,ny) u = np.ones((ny,nx)) ##create a 1xn vector of 1's un = np.ones((ny,nx)) ## ###Assign initial conditions u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2 ##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2 ###Plot Initial Condition fig = plt.figure(figsize=(11,7), dpi=100) ##the figsize parameter can be used to produce different sized images ax = fig.gca(projection='3d') X, Y = np.meshgrid(x,y) surf = ax.plot_surface(X,Y,u[:]) u = np.ones((ny,nx)) u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2 for n in range(nt+1): ##loop across number of time steps un[:] = u[:] for i in range(1, len(u)): for j in range(1, len(u)): u[i,j] = un[i, j] - (c*dt/dx*(un[i,j] - un[i-1,j]))-(c*dt/dy*(un[i,j]-un[i,j-1])) u[0,:] = 1 u[-1,:] = 1 u[:,0] = 1 u[:,-1] = 1 fig = plt.figure(figsize=(11,7), dpi=100) ax = fig.gca(projection='3d') surf2 = ax.plot_surface(X,Y,u[:]) u = np.ones((ny,nx)) u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2 for n in range(nt+1): ##loop across number of time steps un[:] = u[:] u[1:,1:]=un[1:,1:]-(c*dt/dx*(un[1:,1:]-un[0:-1,1:]))-(c*dt/dy*(un[1:,1:]-un[1:,0:-1])) u[0,:] = 1 u[-1,:] = 1 u[:,0] = 1 u[:,-1] = 1 fig = plt.figure(figsize=(11,7), dpi=100) ax = fig.gca(projection='3d') surf2 = ax.plot_surface(X,Y,u[:]) from IPython.display import YouTubeVideo YouTubeVideo('tUg_dE3NXoY') from IPython.core.display import HTML def css_styling(): styles = open("../styles/custom.css", "r").read() return HTML(styles) css_styling()