dsolve and applying boundary conditions using substitutions.Separation of variables is a technique that reduces a partial differential equation to a pair of ordinary differential equations. Solutions of the resulting ODEs are more straightforward than for the PDE. This lecture and computational notebook will help you analyze one such solution to the diffusion equation. To start we will select a problem that is best solved with a Fourier sin series. Making this selection at the start is arbitrary, however the consequences will define the physics of the problem. This is most easily seen in the restrictions on the boundary conditions. For this problem our boundary conditions will be chosen to be compatible with a Fourier sin series.
In [ ]:
%matplotlib inline
import sympy as sp
# You may have to comment out the init_sesison() line if you have not used my session.py file.
sp.init_session(quiet=True)
sp.init_printing()
dum = sp.symbols('dum')
In [ ]:
def initial_condition(x):
return x
def b_m_amplitudes(m, funToProject, center=0.5, width=1.0):
return (2/width)*sp.integrate(funToProject(dum)*sp.sin(2*m*sp.pi*dum/width), (dum,center-width/2,center+width/2))
def b_m_vectorspace_element(m, var, width=1.0):
return sp.sin(2*m*sp.pi*var/width)
terms = 10
width = 2
center = 0
sinSeries = sum([b_m_amplitudes(m, initial_condition, center=center, width=width)*b_m_vectorspace_element(m,x,width=width) for m in range(terms)])
In [ ]:
sp.plot(sinSeries,(x,0,1));
The physical problem of interest is one where an extensive quantity (energy, mass, etc.) evolves diffusivly in some medium. Using mass diffusion as an example, we will identify the concentration of species in our medium at any point and tiime, $c(x,t)$. We must choose a domain (the geometry of the medium) in which this diffusion process will take place. In the solution developed below we assume that we have a long bar of material with a constant cross section and this bar is insulated so that none of our extensive quantity can escape through the ends or through the cylindrical surface. These assumptions permit simplification of the problem to one dimension.
We will assume that our bar of material has a fixed length $L = 1.0$ units long. The medium will contain an initial amount of the diffusing substance whose distribution will be given by a function that represents the initial condition. We will represent our initial condition using the symbol $\phi(x,0)$ where the zero indicates time $t=0$. For this exercise we will choose:
$$ \phi(x,0) = x $$The solution to the problem depends in part on our ability to represent the initial condition as a Fourier series. The choice of $\phi(x)=x$ permits us demonstrate something about the series solution at the boundary. Let us visualize the solution to get a feel for what we're dealing with first and then we will develop the solution ourselves later in this notebook. In the interact function we can change the time, the number of terms, and the diffusion constant ($\alpha$).
solutionAtX that takes a pair of values ($x$, $t$), a diffusion coefficient $\alpha$ and the number of terms to use in the Fourier series.make_plot) that we pass to interact so that we have an interactive plot to examine the solution and a set of parameters.
In [ ]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
def solutionAtX(x, t, alpha, terms):
# This is an already developed solution to the diffusion equation.
# If you can follow along in the lecture - your solution will look
# identical.
return sum([
(2/(np.pi*n)**2)*
(np.sin(n*np.pi*x))*
(-np.pi*n*np.cos(np.pi*n)+np.sin(np.pi*n))*
np.exp(-t*(alpha*n*np.pi)**2)
for n in range(1,terms+1,1)
])
def make_plot(t=0.01, alpha=0.2, terms=20):
plt.figure(figsize=(7,5))
x = np.arange(0, 1, 0.001)
y = solutionAtX(x, t, alpha, terms)
plt.plot(x, y, x, x, 'orange')
plt.fill_between(x, 0, y, alpha=0.2)
plt.ylim(0.0, 1.0)
plt.xlabel('Distance')
plt.ylabel('Concentration')
plt.show()
In [ ]:
interact(make_plot, t=(0.01,2,0.01), alpha=(0.2,0.5,0.01), terms=(1,40,1));
Using the interactive widgets above, explore different combinations of parameters and connect your physical intuition for diffusion with the solution. The orange line is the line $\phi(x,0)=x$ and the blue line is the solution to the diffusion equation. As you advance time try and connect your physical intuition of diffusion with the evolution of the area under the blue curve in the solution. At the boundaries the diffusing species is allowed to escape the medium. Over time, we should see a gradual decrease in the amount of substance (shaded area under the blue curve) in the solution to the equation.
There are many ways to solve a PDE. Separating variables is one such method to reduce a PDE to a pair (or more) of ODEs. We start by assuming that the solution has the form of a time dependent piece and a spatially dependent piece. So we could write:
$$ c(x,t) = X(x) T(t) $$We substitute the assumed form of the solution into Fick's law:
$$ \frac{\partial c(x,t)}{\partial t} = D \frac{\partial^2 c(x,t)}{\partial x^2} $$this gives us:
$$ \frac{\partial X(x) T(t)}{\partial t} = D \frac{\partial^2 X(x) T(t)}{\partial x^2} $$We will differentate the assumed solution of $c(x,t)$ above according to the derivatives in the differential equation. It is usual to let D be $\alpha^2$ to guarantee that D is always positive, we can do this in sympy with a keyword in the symbols function call.
In [ ]:
import sympy as sp
sp.init_session(quiet=True)
sp.init_printing()
In [ ]:
var('X,T',cls=Function)
In [ ]:
var('alpha gamma c', positive=True)
With our substitutions in place:
$$ \frac{\partial X(x) T(t)}{\partial t} = \alpha^2 \frac{\partial^2 X(x) T(t)}{\partial x^2} $$After taking the derivatives we have:
$$ X(x)T'(t) = \alpha^2 X''(x)T(t) $$We can arrange these expressions so that all the functions of "t" are on one side and the functions of "x" are on the other side. Divide both the LHS and RHS of the above expression by $\alpha^2 X(x) T(t)$ to seperate variables on either sides of the equals sign. Note also that this means that the product $\alpha^2 X(x) T(t)$ should never be zero. Add a seperation constant $\gamma$ and integrate.
The rationalization for the seperation constant is that since the LHS and RHS each depend on one variable then the only way the LHS and RHS can be equal is if they are both equal to a constant. Our ODEs are then written as:
In [ ]:
timeODE = sp.Eq(T(t).diff(t,1)/(alpha**2*T(t)),-gamma**2)
timeODE
In [ ]:
spaceODE = sp.Eq(X(x).diff(x,2)/X(x),-gamma**2)
spaceODE
Invoke dsolve()
In [ ]:
sp.dsolve(timeODE,T(t))
This is not the form you are likely to find in a textbook. It is the correct solution, but it may not be the best solution for us. We can examine all the possible forms of the solution so that we can inspect the list and pick the one we like. We create a list of hints and then use a list comprehension to solve the equation for each hint. The notebook displays all the solutions so that we can find the one that suits our needs.
In [ ]:
list_of_hints = sp.classify_ode(timeODE)
list_of_solutions = [sp.dsolve(timeODE,T(t), hint=hint) for hint in list_of_hints]
list_of_solutions
A bit of matching up with the list of hints and we find that the one we are interested in is 1st_linear:
In [ ]:
timePiece = sp.dsolve(timeODE,T(t), hint='1st_linear')
timePiece
In [ ]:
spacePiece = sp.dsolve(spaceODE,X(x))
spacePiece
The solution (with constants) is the product of the two solutions above:
$$ c(x,t) = X(x)T(t) = \Big( C_{1} \sin{\left (\gamma x \right )} + C_{2} \cos{\left (\gamma x \right )} \Big)\Big ( C_{3} e^{- \alpha^{2} \gamma^{2} t} \Big ) $$We will need to "absorb" the time solution constant ($C_3$) into the space solution constants. There are an infinity of solutions depending on our choices of $C_1$ and $C_2$. These will ultimately be determined by the boundary conditions.
I would also note here the following: the CAS (SymPy in this case) did all the work for us above. The code below is really to just keep ourselves organized. I don't intend to do much with generalSolution except to make substitutions to find the boundary conditions. Working with a CAS is equal parts setup, organization, and deciphering the results.
Clean up the solution for use below. Eliminate one of the constants as the product of C1 and C2 is also a constant. Sympy does not automatically create the symbols corresponding to the constants for us. So we need to create those before proceeding:
In [ ]:
var('C1 C2')
In [ ]:
generalSolution = Eq(c(x,t),(timePiece.rhs).subs([(C1,1)])*spacePiece.rhs)
generalSolution
The following conditions apply to this problem:
$$ c(0,t) = 0 $$$$ c(1,t) = 0 $$$$ c(x,0) = x $$These conditions apply for the left-hand boundary, the right-hand boundary and the initial condition (the time boundary). The domain of the experiment is $\{\,x\, \mid 0 < x < 1\,\}$
In [ ]:
generalSolution
Start by substituting the boundary condition $c(0,t) = 0$:
In [ ]:
generalSolution.subs({x:0})
The only way this can be true is if $C_2$ is zero. Make that substition and insert the second boundary condition where $c(1,t)$ must also be zero:
In [ ]:
generalSolution.subs({x:1, C2:0})
So - wait, what? Can can we let $C1=0$ too?
That IS a solution. It just is not the one we're looking for, so we need to examine how the third factor (the sin term) can help us.
What about the exponential? At the boundary there is no dependence of $x$ in this factor however we know that the parameters in the exponential will not cause that factor to go to zero. It starts at one and just gets smaller - so looking at the second factor will not help us much.
This means that we have to let $\sin(\gamma)=0$. We know this can happen with the right choice of $\gamma$. What choices can we use?
What if:
$$\gamma = \pm \pi, \pm 2\pi, \pm 3\pi, ... $$and in general:
$$\gamma = \pm n\pi$$Let us go back an examine the initial condition where $t=0$.
In [ ]:
generalSolution.subs({t:0, C2:0, gamma:n*sp.pi})
Now this statement tells us that we need to find a constant $C_1$ such that $c(x,0)=x$. You may be wondering how we can do that - but you already know how to do that.
We identify the diffusion PDE (Fick's Law) as a homogeneous PDE. Kreyszig defines homogeneous as an equation in $u$ where "each of its terms contains either $u$ or one of its partial derivatives". Because the equation is a homogeneous PDE we can invoke the principle of superposition that combinations of solutions to the PDE are also solutions to the PDE.
So it is therefore possible to write (without proof) that the initial condition:
$$c(x,0) = \sum_{n=1}^N C_n sin(n\pi x) = \phi(x) = x$$and the general solution to the equation is:
$$c(x,t) = \sum_{n=1}^N C_n sin(n\pi x) e^{- \alpha^2 n^2 \pi^2 t}$$The summation is the Fourier sine series - so as long as you can compute the coefficients $C_n$ then you can write down a summation that will be an approximate solution to the problem. All that remains is to find the Fourier coefficient for the function $\phi(x)$ and substitute those into the general solution. Go back and look at the plotted solution and see if you can write your own function derived from the general solution and a Fourier series.
Your homework assignment is to set up a different initial condition and solve the problem again.
From Farlow's book entitled "Partial Differential Equations for Scientists and Engineers" he poses the following two problems:
sin series you will also be requiring the values of the function to be zero at the boundaries. Hint: The domain of the series may be larger than the domain of the physical problem.Note that there is a physical impossibility in this problem. Can you identify it and discuss as part of your solution?