Open Educational Resources

Partial Differential Equations: Semi-Discretization

The following code is compiled in the Jupyter Notebook.

Example 13.6a

View Python Code

import numpy as np
from numpy import pi,sin,exp,diag
import matplotlib.pyplot as plt
%matplotlib notebook
# Physical parameters
alpha = 1
L = 1
# Mesh
N = 21
dx = 1/(N-1)
x = np.linspace(0,L,N)
# Time advancement
dt = 0.00125
tmax = 2
t = np.arange(0,tmax+dt/10,dt)
nmax = t.size
T = np.array([[0.0 for i in range(N)] for j in range(nmax+1)])
# initial condition
T[0] = sin(pi*x)
# inhomogeneous part
fvec = lambda x,t: (pi**2-1)*exp(-t)*sin(pi*x)
# Matrix
twos = np.array([-2 for i in range(N-2)])
ones = np.array([1 for i in range(N-3)])
A = diag(twos)+diag(ones,k=1)+diag(ones,k=-1)
# Explicit Euler
for n in range(nmax):
    dTdt = (alpha/dx**2)*A@T[n][1:N-1]+fvec(x[1:N-1],t[n])
    T[n+1][1:N-1] = T[n][1:N-1] + dt*dTdt
fig,ax = plt.subplots(1,1)
ax.grid(True, which='both')
plt.title("Explicit Euler solution of Heat Equation with dt = " + str(dt))
plt.xlabel("x")
plt.ylabel("u")
ax.plot(x,T[0])
for n in range(nmax):
    ax.lines[0].set_ydata(T[n])
    fig.canvas.draw()
fig,ax = plt.subplots(1,1)
ax.grid(True, which='both')
plt.title("Explicit Euler solution of Heat Equation with dt = " + str(dt))
plt.xlabel("x")
plt.ylabel("u")
for n in range(0,nmax,400):
    ax.plot(x,T[n], label = "t =" + str(t[n]))
    ax.legend()
    fig.canvas.draw()

Example 13.6b

View Python Code

import numpy as np
from numpy import pi,sin,exp,diag
import matplotlib.pyplot as plt
%matplotlib notebook
# Physical parameters
alpha = 1
L = 1

# Mesh
N = 21
dx = 1/(N-1)
x = np.linspace(0,L,N)

# Time advancement
dt = 0.0015
tmax = 2
t = np.arange(0,tmax+dt/10,dt)
nmax = t.size
T = np.array([[0.0 for i in range(N)] for j in range(nmax+1)])
T[0] = sin(pi*x)
fvec = lambda x,t: (pi**2-1)*exp(-t)*sin(pi*x)
# Matrix
twos = np.array([-2 for i in range(N-2)])
ones = np.array([1 for i in range(N-3)])
A = diag(twos)+diag(ones,k=1)+diag(ones,k=-1)
for n in range(nmax):
    dTdt = (alpha/dx**2)*A@T[n][1:N-1]+fvec(x[1:N-1],t[n])
    T[n+1][1:N-1] = T[n][1:N-1] + dt*dTdt
fig,ax = plt.subplots(1,1)
ax.grid(True, which='both')
plt.title("Explicit Euler solution of Heat Equation with dt = " + str(dt))
plt.xlabel("x")
plt.ylabel("u")
ax.plot(x,T[0])
for n in range(nmax):
    ax.lines[0].set_ydata(T[n])
    fig.canvas.draw()
fig,ax = plt.subplots(1,1)
ax.grid(True, which='both')
plt.title("Explicit Euler solution of Heat Equation with dt = " + str(dt))
plt.xlabel("x")
plt.ylabel("u")
tplots = [0,100,110,114,118]
for n in tplots:
    ax.plot(x,T[n], label = "t =" + str(t[n]))
    ax.legend()
    fig.canvas.draw()

Example 13.7a

View Python Code

import numpy as np
from numpy import pi,sin,exp,diag
import matplotlib.pyplot as plt
%matplotlib notebook
# Physical parameters
c = 1
L = 0.7

# Mesh
dx = 0.01
N = int(L/dx +1)
x = np.arange(0,L+dx/10,dx)
x
# Time advancement
dt = 0.01
tmax = 0.6
t = np.arange(0,tmax+dt/10,dt)
nmax = t.size
u = np.array([[0.0 for i in range(N)] for j in range(nmax+1)])
# initial condition
u[0] = exp(-200*(x-0.25)**2)
# Define Matrix
zeros = np.array([0 for i in range(N-1)])
ones = np.array([1 for i in range(N-2)])
A = -1*(c/2/dx)*(diag(zeros)+diag(ones,k=1)-diag(ones,k=-1))
A[N-2,N-3] = c/dx
A[N-2,N-2] = -c/dx
# Explicit Euler
for n in range(nmax):
    dudt = A@u[n][1:N]
    u[n+1][1:N] = u[n][1:N] + dt*dudt
fig,ax = plt.subplots(1,1)
ax.grid(True, which='both')
plt.title("Explicit Euler solution of Convection Equation with dt = " + str(dt))
plt.xlabel("x")
plt.ylabel("u")
ax.plot(x,u[0])
for n in range(nmax):
    ax.lines[0].set_ydata(u[n])
    fig.canvas.draw()
fig,ax = plt.subplots(1,1)
ax.grid(True, which='both')
plt.title("Explicit Euler solution of Convection Equation with dt = " + str(dt))
plt.xlabel("x")
plt.ylabel("u")
tplots = [0,12,25,38]
for n in tplots:
    ax.plot(x,u[n], label = "t =" + str(t[n]))
    ax.legend()
    fig.canvas.draw()

Example 13.7b

View Python Code

import numpy as np
from numpy import pi,sin,exp,diag
import matplotlib.pyplot as plt
%matplotlib notebook
# Physical parameters
c = 1
L = 1

# Mesh
dx = 0.01
N = int(L/dx +1)
x = np.arange(0,L+dx/10,dx)
x
# Time advancement
dt = 0.01
tmax = 1.1
t = np.arange(0,tmax+dt/10,dt)
nmax = t.size
u = np.array([[0.0 for i in range(N)] for j in range(nmax+1)])
# initial condition
u[0] = exp(-200*(x-0.25)**2)
# Define Matrix
zeros = np.array([0 for i in range(N-1)])
ones = np.array([1 for i in range(N-2)])
A = -1*(c/2/dx)*(diag(zeros)+diag(ones,k=1)-diag(ones,k=-1))
A[N-2,N-3] = c/dx
A[N-2,N-2] = -c/dx
for n in range(nmax):
    k1 = A@u[n,1:N]
    k2 = A@u[n,1:N] + dt/2.*A@k1
    k3 = A@u[n,1:N] + dt/2.*A@k2
    k4 = A@u[n,1:N] + dt*A@k3
    u[n+1,1:N] = u[n,1:N] + dt*(k1+2*k2+2*k3+k4)/6  

fig,ax = plt.subplots(1,1)
ax.grid(True, which='both')
plt.title("RK4 solution of Convection Equation with dt = " + str(dt))
plt.xlabel("x")
plt.ylabel("u")
ax.plot(x,u[0])
for n in range(nmax):
    ax.lines[0].set_ydata(u[n])
    fig.canvas.draw()
fig,ax = plt.subplots(1,1)
ax.grid(True, which='both')
plt.title("RK4 solution of Convection Equation with dt = " + str(dt))
plt.xlabel("x")
plt.ylabel("u")
tplots = [0,25,50,75]
for n in tplots:
    ax.plot(x,u[n], label = "t =" + str(t[n]))
    ax.legend()
    fig.canvas.draw()

Leave a Reply

Your email address will not be published. Required fields are marked *