Question: 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
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")
Step by Step Solution
There are 3 Steps involved in it
Get step-by-step solutions from verified subject matter experts
