Boundary Value Problems, Partial Differential Equations, and how to plot our problems in 3D.
That may sound like a lot, however, it is very similar to what you have already learned. If a partiticular section is difficult to understand, try looking at the code line by line.
Our standard matplotlib cannot handle 3D graphing. Just as before we will import
the tools we need.
from mpl_toolkits.mplot3d import Axes3D
We set our Z axes equal to certain data, perhaps a function $z = f(x,y)$
Z = data
We can use numpy.meshgrid
to create our meshgrid
X, Y = np.meshgrid(np.arange(3), np.arange(3))
We create a figure
fig = plt.figure()
Add our plot to the figure
ax = fig.add_subplot(1,1,1, projection='3d')
Choose the arrangement
ax.plot_surface(X, Y, Z)
And that's it!
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data = [[1, 2, 3], [1, 2, 3], [1, 2, 3]]
T_array = [1, 2, 3]
R_array = [3, 4,5]
data = np.array(data)
Z = data
X, Y = np.meshgrid(np.arange(3), np.arange(3))
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
ax.plot_surface(X, Y, Z)
plt.show()
The toolkit does not strictly prevent you from swapping the axes (or their signs)
ax.plot_surface(Y, Z, X)
ax.plot_surface(-X, -Z, Y)
However, if you are inconsistent, you risk confusing data and graphs. Be especially careful when plotting multiple things onto the same figure.
Here two planes will appear on the same subplot, and the X of one is the Y of another. You are not prevented from doing this.
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
ax.plot_surface(Y, Z, X)
ax.plot_surface(X, Z, Y)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data = [[1, 2, 3], [1, 2, 3], [1, 2, 3]]
T_array = [1, 2, 3]
R_array = [3, 4,5]
data = np.array(data)
Z = data
X, Y = np.meshgrid(np.arange(3), np.arange(3))
fig1 = plt.figure()
ax1 = fig1.add_subplot(1,1,1, projection='3d')
ax1.plot_surface(Y, Z, X)
fig2 = plt.figure()
ax2 = fig2.add_subplot(1,1,1, projection='3d')
ax2.plot_surface(-X, -Z, Y)
plt.show()
Optionally, just as in 2D, we can add labels:
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
Unfortunately, the $z$ label does not automatically fit within the plot.
Furthermore, the labels assume you have chosen the default:
ax.plot_surface(X, Y, Z)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data = [[1, 2, 3], [1, 2, 3], [1, 2, 3]]
T_array = [1, 2, 3]
R_array = [3, 4,5]
data = np.array(data)
Z = data
X, Y = np.meshgrid(np.arange(3), np.arange(3))
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
ax.plot_surface(X, Y, Z)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
plt.show()
We have a function, $z = f(x,y) = x^2 - y^2$, where $-L < x < L$ and $-L < y < L$, provided $L = 10$, give the 3D plot
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# https://en.wikipedia.org/wiki/Saddle_point
L = 10.
x = np.linspace(-L, L, 100)
y = np.linspace(-L, L, 100)
X, Y = np.meshgrid(x, y)
Z = (X**2 - Y**2)
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
ax.plot_surface(X, Y, Z);
# plt.show();
Of course the number of intervals affects the resolution and accuracy:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# https://en.wikipedia.org/wiki/Saddle_point
L = 10.
x = np.linspace(-L, L, 8)
y = np.linspace(-L, L, 8)
X, Y = np.meshgrid(x, y)
Z = (X**2 - Y**2)
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
ax.plot_surface(X, Y, Z);
# plt.show();
Before we step into the territory of code, let’s find out what Partial Differential Equation
(PDE) means, and how it differs from Ordinary Differential Equations
(ODE).
Wikipedia explains this well, whereas an ODE “is a differential equation containing one or more functions of one independent variable and its derivatives”, a PDE “may be with respect to more than one independent variable”.
In essence, we are doing what we did with ODEs, but simply with more variables.
As you may very well suspect, we solve PDEs in a similar way to ODEs.
There is a specific case of differential equations, where we have an Initial Value Problem (IVP); in such a problem we may have the general solution, but we still need a particular solution in order to solve our problem. It is with the Initial Condition(s) (IC) that we do so.
In layman’s terms, an IVP is a problem where we have a sense of what the function is (either by having calculated it ourselves or by having been provided the function), and our job is to specify the terms of the function, so it can satisfy some conditions.
Now let's start coding a problem involving PDE and ICs. This particular problem is called a Heat Equation
, Diffusion Equation
, or Heat Diffusion Equation
, it describes the diffusion of heat over an object as time progresses.
We define a function, in this case:
def calc_heat(nu):
nx = 41
dx = 2 / (nx - 1)
nt = 20 #the number of timesteps we want to calculate
# nu = 0.3 #the value of viscosity
sigma = .2 #sigma is a parameter, we'll learn more about it later
dt = dx**2 #dt is defined using sigma ... more later!
u = np.ones(nx) #a numpy array with nx elements all equal to 1.
u[int(.5 / dx):int(1 / dx + 1)] = 2 #setting u = 2 between 0.5 and 1 as per our I.C.s
plt.plot(np.linspace(0, 2, nx), u,'k--');
Our heat diffusion has a range of possible forms it can take on depending on the viscosity
un = np.ones(nx) #our placeholder array, un, to advance the solution in time
for n in range(nt): #iterate through time
un = u.copy() #copy the existing values of u into un
for i in range(1, nx - 1):
u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2 * un[i] + un[i-1])
u[0] = 1
u[-1] = 1
plt.plot(np.linspace(0, 2, nx), u,'r--');
import numpy as np #loading our favorite library
import matplotlib.pyplot as plt #and the useful plotting library
plt.style.use('ggplot')
%matplotlib inline
from ipywidgets import interact
from IPython.display import clear_output, display, HTML
def calc_heat(nu):
nx = 41
dx = 2 / (nx - 1)
nt = 20 #the number of timesteps we want to calculate
# nu = 0.3 #the value of viscosity
sigma = .2 #sigma is a parameter, we'll learn more about it later
# dt = sigma * dx**2 / nu #dt is defined using sigma ... more later!
dt = dx**2 #dt is defined using sigma ... more later!
u = np.ones(nx) #a numpy array with nx elements all equal to 1.
u[int(.5 / dx):int(1 / dx + 1)] = 2 #setting u = 2 between 0.5 and 1 as per our I.C.s
plt.plot(np.linspace(0, 2, nx), u,'k--');
un = np.ones(nx) #our placeholder array, un, to advance the solution in time
for n in range(nt): #iterate through time
un = u.copy() #copy the existing values of u into un
for i in range(1, nx - 1):
u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2 * un[i] + un[i-1])
u[0] = 1
u[-1] = 1
plt.plot(np.linspace(0, 2, nx), u,'r--');
plt.xlabel('x')
plt.ylabel('T')
plt.xlim(0,2)
plt.ylim(1,2)
plt.show()
return None
interact(calc_heat, nu=(0., .5, 0.02));
Here is our heat diffusion specifically when $nu = 0.3$
import numpy #loading our favorite library
from matplotlib import pyplot #and the useful plotting library
%matplotlib inline
nx = 41
dx = 2 / (nx - 1)
nt = 20 #the number of timesteps we want to calculate
nu = 0.3 #the value of viscosity
sigma = .2 #sigma is a parameter, we'll learn more about it later
dt = sigma * dx**2 / nu #dt is defined using sigma ... more later!
u = numpy.ones(nx) #a numpy array with nx elements all equal to 1.
u[int(.5 / dx):int(1 / dx + 1)] = 2 #setting u = 2 between 0.5 and 1 as per our I.C.s
un = numpy.ones(nx) #our placeholder array, un, to advance the solution in time
for n in range(nt): #iterate through time
un = u.copy() #copy the existing values of u into un
for i in range(1, nx - 1):
u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2 * un[i] + un[i-1])
pyplot.plot(numpy.linspace(0, 2, nx), u);
np.full_like(x, 673)
Let's continue the previous example in 3D.
Starting with the setup of our tools:
import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D ##library for 3d projection plots
%matplotlib inline
Now we declare our variables:
$nt$ our timesteps, $nu$ our viscosity, $\sigma$, $dx$, and $dt$
Using $nx$ and $ny$, the end of our ranges, we define two linspaces:
x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 2, ny)
And also our initial conditions (ICs).
###variable declarations
nx = 31
ny = 31
nt = 17
nu = .01
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .25
dt = sigma * dx * dy / nu
x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 2, ny)
u = numpy.ones((ny, nx)) # create a 1xn vector of 1's
un = numpy.ones((ny, nx))
###Assign initial conditions
# set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2
With our linspaces we define a meshgrid:
x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 2, ny)
X, Y = numpy.meshgrid(x, y)
And we graph as before,
starting with a figure
fig = pyplot.figure()
Add our plot to the figure
ax = fig.gca(projection='3d')
Choose the arrangement
surf = ax.plot_surface(X, Y, u, rstride=1, cstride=1, cmap=cm.viridis, linewidth=0, antialiased=False)
ax.set_xlim(0, 2)
ax.set_ylim(0, 2)
ax.set_zlim(1, 2.5)
fig = pyplot.figure()
ax = fig.gca(projection='3d')
X, Y = numpy.meshgrid(x, y)
surf = ax.plot_surface(X, Y, u, rstride=1, cstride=1, cmap=cm.viridis,
linewidth=0, antialiased=False)
ax.set_xlim(0, 2)
ax.set_ylim(0, 2)
ax.set_zlim(1, 2.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$');
Just like in the 2D example, we are going to simulate heat diffusion, there is a lot here, so let's go through this step by step.
We define our diffuse function, which takes in the timestep as a parameter:
def diffuse(nt):
We define a 2D array $u$ which is defined by the derivatives of $y$ and $x$
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2
We iterate through time, and copy the existing values of $u$ into $un$, just like before. However, this time we have an extra dimension, so we iterate over every $i$ and $j$ within $u$.
for n in range(nt + 1):
un = u.copy()
for i in range(1, u.shape[0]-1):
for j in range(1, u.shape[1] - 1):
u[i,j] = un[i, j] + nu * dt / dx**2 * \
(un[i+1, j] - 2*(un[i,j]) + un[i-1,j]) + \
nu * dt / dy**2 * \
(un[i,j+1] - 2*un[i, j] + un[i, j-1])
Here we are setting the outer square of our 2D array to 1, convince yourself that this is correct:
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1
And now we plot the 3D graph.
-Starting with a figure
-Add our plot to the figure
-Choose the arrangement
fig = pyplot.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis,
linewidth=0, antialiased=True)
ax.set_zlim(1, 2.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$');
###Run through nt timesteps
def diffuse(nt):
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2
# un+1i,j=uni,j+νΔtΔx2(uni+1,j−2uni,j+uni−1,j)+νΔtΔy2(uni,j+1−2uni,j+uni,j−1)
# u[1:-1, 1:-1] = (un[1:-1,1:-1] +
# nu * dt / dx**2 *
# (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
# nu * dt / dy**2 *
# (un[2:,1: -1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))
for n in range(nt + 1):
un = u.copy()
for i in range(1, u.shape[0]-1):
for j in range(1, u.shape[1] - 1):
u[i,j] = un[i, j] + nu * dt / dx**2 * \
(un[i+1, j] - 2*(un[i,j]) + un[i-1,j]) + \
nu * dt / dy**2 * \
(un[i,j+1] - 2*un[i, j] + un[i, j-1])
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1
fig = pyplot.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis,
linewidth=0, antialiased=True)
ax.set_zlim(1, 2.5)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$');
# http://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/09_Step_7.ipynb
Now we plug in 10 timesteps into the diffuse function we defined before:
diffuse(10)
Let’s take another moment before we start coding again. What is a Boundary Value Problem (BVP)? Well, it is a type of problem where we must satisfy certain boundary conditions. If this sounds familiar, it should be. Our old friend Wikipedia comes to save the day again, “a [BVP] has conditions specified at the extremes (‘boundaries’) of the independent variable in the equation whereas an [IVP] has all of the conditions specified at the same value of the independent variable”.
Imagine we have some problem involving the height of a skydiver
an IVP would define a function $h(t)$,
perhaps involving $h'(t) = k \cdot h(t)$ and $h(0) = C$, where k and C are specified
Here the we only ever use the value $t_i = 0$ in $h(0)$
How about a BVP? Well, here there has been defined two boundaries
, such as the initial height $h(0) = h_1$ and the final height $h(t_f) = h_2$.
We would need the differential equation as well. However, you can see here we use $t_i$ and $t_f$. These were the aforementioned "extremes", before the skydiver jumps (the highest possible height) and after they land (the lowest possible height).
For our heat equations
, these boundaries could be end points on an object. There are multitude of problems that need boundary conditions in order to solve.
We have a few new tools here as well:
from scipy.optimize import least_squares
from scipy.integrate import ode
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
# import snips as snp
# snp.prettyplot(matplotlib)
from scipy.optimize import least_squares
from scipy.integrate import ode
from scipy.integrate import odeint
%matplotlib inline
The method we will be using is solve by shooting.
We will test our ODE function, "odefun", and see if overshoots our initial boundary condition or undershoots it.
Provided that we know a limit $u_1(0)$, our miss-shot can give us an $u_2(0)$ using fsolve.
In this case we will guess $u_2(0) = 1$
Our odefun is:
def odefun(U, y):
u1, u2 = U
mu = 1
Pdrop = -100
du1dy = u2
du2dy = 1.0 / mu * Pdrop
return [du1dy, du2dy]
# We need u_1(0) and u_2(0), but we only have u_1(0). We need to guess a value for u_2(0)
# and see if the solution goes through the u_2(d)=0 boundary value.
d = 0.1 # plate thickness
def odefun(U, y):
u1, u2 = U
mu = 1
Pdrop = -100
du1dy = u2
du2dy = 1.0 / mu * Pdrop
return [du1dy, du2dy]
u1_0 = 0 # known
u2_0 = 1 # guessed
dspan = np.linspace(0, d)
U = odeint(odefun, [u1_0, u2_0], dspan)
plt.plot(dspan, U[:,0])
plt.plot([d],[0], 'ro')
plt.xlabel('d')
plt.ylabel('$u_1$');
# plt.savefig('images/bvp-shooting-1.png')
And so we see here that we have undershot, we want to overshoot for our range, so we try again. This time with $u_2(0) = 7$.
You can interact with the $u_2(0)$ here to see how it affects our odefun
from ipywidgets import interact, interactive, fixed, interact_manual
def calc_poiss(u2_0):
u1_0 = 0 # known
# u2_0 = 7 # guessed
dspan = np.linspace(0, d)
U = odeint(odefun, [u1_0, u2_0], dspan)
plt.plot(dspan, U[:,0])
plt.plot([d],[0], 'ro')
plt.xlabel('d')
plt.ylabel('$u_1$')
interact(calc_poiss, u2_0=(0., 7, .2));
From the previous slider, you can see that we satisfy our boundary condition when $u_2(0) \simeq 5$, but we want to know how to calculate this number, and we want something more precise.
as mentionned earlier, since we know $1 < u_2(0) < 7$ we can use fsolve to figure this out for us. Luckily that comes with scipy
from scipy.optimize import fsolve
And we have our known initial condition, and the size of d spans from 0 to 0.1
u1_0 = 0 # known
dspan = np.linspace(0, d)
Now we need to give fsolve an function to solve. This is our objective function.
def objective(u2_0):
U = odeint(odefun, [u1_0, u2_0], dspan)
u1 = U[:,0]
return u1[-1]
Now we set $u_2(0)$ equal to what we determined it should be.
u2_0, = fsolve(objective, 1.0)
And using odeint as before, we get our Numerical solution:
U = odeint(odefun, [u1_0, u2_0], dspan)
plt.plot(dspan, U[:,0], label='Numerical solution')
plt.plot([d],[0], 'ro')
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
d = 0.1 # plate thickness
Pdrop = -100
mu = 1
def odefun(U, y):
u1, u2 = U
du1dy = u2
du2dy = 1.0 / mu * Pdrop
return [du1dy, du2dy]
u1_0 = 0 # known
dspan = np.linspace(0, d)
def objective(u2_0):
U = odeint(odefun, [u1_0, u2_0], dspan)
u1 = U[:,0]
return u1[-1]
u2_0, = fsolve(objective, 1.0)
# now solve with optimal u2_0
U = odeint(odefun, [u1_0, u2_0], dspan)
plt.plot(dspan, U[:,0], label='Numerical solution')
plt.plot([d],[0], 'ro')
# plot an analytical solution
u = -(Pdrop) * d**2 / 2 / mu * (dspan / d - (dspan / d)**2)
plt.plot(dspan, u, 'r--', label='Analytical solution')
print(u2_0)
plt.xlabel('d')
plt.ylabel('$u_1$')
plt.legend(loc='best');