Homework SoursePython plots I have this script for a project However when I
Python plots:
I have this script for a project. However, when I run the script the plots are not showing. How do I get the plots to display?
from math import pi
import pylab
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from numpy import linspace, exp, pi, 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()
plt.xlabel(\"x\")
plt.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)
plt.xlabel(\"Time\")
plt.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)
plt.xlabel(\"Expansion Terms, N\")
plt.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, \'.\')
plt.xlabel(\"n\")
plt.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)
plt.xlabel(\"Time\")
plt.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))
plt.xlabel(\"x\")
plt.ylabel(\"u\")
plt.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()
plt.xlabel(\"x\")
plt.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()
plt.xlabel(\"x\")
plt.ylabel(\"time\")
Solution
import the following .............
**************************
import matplotlib.pyplot as plt
import numpy as np
****************************
example:
import matplotlib.pyplot as plt
import numpy as np
s = np.arange(0.0, 2.0, 0.01)
t = 1 + np.sin(2*np.pi*s)
plt.plot(s, t)
plt.xlabel(\'time (t)\')
plt.ylabel(\'voltage (mV)\')
plt.title(\'About as simple as it gets, folks\')
plt.grid(True)
plt.savefig(\"test.png\")
plt.show()
example:
import numpy as np
import matplotlib.pyplot as plt
def f(t):
return np.exp(-t) * np.cos(2*np.pi*t)
t1 = np.arange(0.0, 5.0, 0.1)
t2 = np.arange(0.0, 5.0, 0.02)
plt.figure(1)
plt.subplot(211)
plt.plot(t1, f(t1), \'bo\', t2, f(t2), \'k\')
plt.subplot(212)
plt.plot(t2, np.cos(2*np.pi*t2), \'r--\')
plt.show()
code:
import matplotlib.pyplot as plt
import numpy as np
# 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()
plt.xlabel(\"x\")
plt.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)
plt.xlabel(\"Time\")
plt.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)
plt.xlabel(\"Expansion Terms, N\")
plt.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, \'.\')
plt.xlabel(\"n\")
plt.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)
plt.xlabel(\"Time\")
plt.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))
plt.xlabel(\"x\")
plt.ylabel(\"u\")
plt.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()
plt.xlabel(\"x\")
plt.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()
plt.xlabel(\"x\")
plt.ylabel(\"time\")
