Question: Python coding. I'm having trouble running a code a professor gave for a project. When I copy and paste this and run exactly, I am
Python coding. I'm having trouble running a code a professor gave for a project. When I copy and paste this and run exactly, I am getting a bunch of errors. Where do the proper indentions need to go to be able to run this code?
import pylab from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
# Discretization N = 100 nX = 200 nT = 100
# Physical Constants alpha = 1
# T(t) def T(L,c,t): return exp(-L*L * c * t)
# A_n for sine series def A(m): if(m==2): return 0.5
else: return -4*sin(0.5*pi*m)/(pi*m**2 - 4*pi)
# X(x) def X(m,x): return sin(m*pi*x) def dX(m,x): return m*pi*cos(m*pi*x) def ddX(m,x): return -(m*pi)**2*sin(m*pi*x)
# u = X(x)T(t) def U(alpha,x,t,N): Ys = 0 for n in range(1,N+1): Ys = Ys + T(n*pi,alpha,t)*A(n)*X(n,x) return Ys
# du/dx = dX(x)/dx T(t) def dU(alpha,x,t,N): Ys = 0 for n in range(1,N+1): Ys = Ys + T(n*pi,alpha,t)*A(n)*dX(n,x) return Ys
# du/dt = ddX(x)/ddx T(t) def ddU(alpha,x,t,N): Ys = 0 for n in range(1,N+1):
Ys = Ys + T(n*pi,alpha,t)*A(n)*ddX(n,x) return Ys
# Create an x-t grid (appears as a matrix)
[Xm,Tm] = np.mgrid[0:1:nX*1j, 0:0.2:nT*1j]
# Solve U at each x-t point Um = U(alpha,Xm,Tm,N)
def plotXTUsurf(): # Plot the x-t-u surface and a contour above the surface fig = plt.figure() ax = fig.gca(projection='3d') ax.plot_surface(Xm, Tm, Um, rstride=2, cstride=2, alpha=0.9, facecolors=cm.jet(Um)) cset = ax.contour(Xm, Tm, Um, zdir='z', offset=1, cmap=cm.coolwarm) plt.show() xlabel("x") ylabel("time")
# Plot the temperature response at x=0.6 def plotTempX0p6(): t = linspace(0,1,1000) Ut = U(alpha,0.6,t,50)
fig2 = plt.figure() ax2 = fig2.gca() ax2.plot(t, Ut) xlabel("Time") ylabel("u")
print "Maximum Value at x =0.6 is " + str(max(Ut)) print "Located at " + str(t[argmax(Ut)]) # L squared error versus N def L2(N):
Error = 0 for n in range(1,N+1): Error += -1.0/2.0*A(n)**2 Error += 0.25 return Error
def PlotL2Error(): n = zeros(N) L2E = zeros(N) for i in range(1,N): n[i] = i+1
L2E[i] = L2(i+1) fig3 = plt.figure() ax3 = fig3.gca() ax3.semilogy(L2E) xlabel("Expansion Terms, N") ylabel("$L_2$ Error")
# Coefficients def plotCoeffs(): ACoeffs = zeros(N) for i in range(1,N): ACoeffs[i]=A(i) fig4 = plt.figure() ax4 = fig4.gca() ax4.semilogy(ACoeffs,'.') xlabel("n") ylabel("$A_n$")
def FiniteDiffCheck(t,x,h,N): dt = ( U(alpha,x,t+h,N) - U(alpha,x,t-h,N)) /(2*h) dxx = ( U(alpha,x+h,t,N) - 2.0* U(alpha,x,t,N) + U(alpha,x-h,t,N)) /(h**2) return dt-alpha*dxx print(FiniteDiffCheck(0.25,0.25,0.001,100))
# Calculate the energy the hard way # Create array of du/dx on both boundaries def TotalEnergy(): TimeEffectivelyZero = 4 t = linspace(0,TimeEffectivelyZero,5000) dUtLeft = dU(alpha,0,t,50) dUtRight = dU(alpha,1.0,t,50) fig5 = plt.figure() ax5 = fig5.gca() ax5.plot(t, dUtLeft) ax5.plot(t, dUtRight) xlabel("Time") ylabel("u")
# Create a rough trapezoidal rule integration routine def TrapIntegrate(array,a,b): n = size(array) I=0
for i in range(0,n-1): I = I + 0.5*(array[i]+array[i+1]) return I*(b-a)/(n-1) #Energy via flux print("Energy Transport via Fluxes " + str(TrapIntegrate(dUtLeft-dUtRight,0,TimeEffectivelyZero)))
# Check Initial Condition def CheckIC(): x = linspace(0,1,300) fig6 = plt.figure()
ax6 = fig6.gca() for i in [1,2,4,8,16,32,64,128,256,512]: U0 = U(alpha,x,0,i) ax6.plot(x, U0, label="N="+str(i)) xlabel("x") ylabel("u") legend(loc='upper right')
# Plot the x-t-u surface and a contour above the surface def FastestChange(): dT = ddU(alpha,Xm,Tm,N) fig7 = plt.figure()
ax7 = fig7.gca(projection='3d') ax7.plot_surface(Xm, Tm, dT, rstride=2, cstride=2, alpha=0.9, facecolors=cm.jet(Um)) cset = ax7.contour(Xm, Tm, dT, zdir='z', offset=1, cmap=cm.coolwarm) plt.show() xlabel("x") ylabel("time")
# Plot the x-t-u surface and a contour above the surface def ZeroFlux(): dT = dU(alpha,Xm,Tm,N) fig7 = plt.figure()
ax7 = fig7.gca() V = linspace(-1,1,10) cset = ax7.contourf(Xm, Tm, Um,V, cmap=cm.jet) V2 = [0] cset = ax7.contour(Xm, Tm, dT, V2, cmap=cm.binary_r) plt.show() xlabel("x") ylabel("time")
Step by Step Solution
There are 3 Steps involved in it
Get step-by-step solutions from verified subject matter experts
