Homework SoursePython coding Im having trouble running a code a professor g
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\")
Solution
# here is code with no error.
# its output is -4.55316801967e-06
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from math import pi,exp,sin
# 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\")
