# Matplotlib compatibility patch for Pyodide
import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
matplotlib.RcParams._get = dict.get
1.2. Diffusion Equation#
This section builds on Chapter Numerical Modelling and discusses the numerical solution of the 1D diffusion equation.
Note
The learning objectives of this section are:
describe the role of diffusion and identify the diffusion equation
formulate and discretize the diffusion equation using the finite difference method
compute the truncation error for a given numerical scheme
explain the notions of consistency and convergence
explain the notion of stability
derive the stability condition for a conditionally stable scheme
Mathematical model#
The (initial) boundary value problem (IBVP) we wish to solve is the following
The diffusion equation is also known as the heat equation. This equation describes the flow of heat in a rod with a finite length, made out of some heat-conducting material, subject to some boundary conditions at each end. Specifically, this heat spreads out spatially as time increases. This phenomenon is called diffusion. The solution to the heat equation is the temperature distribution \(T(x,t)\) at any time \(t\) and vary with \(x\) along the rod.
We assume that \(T(x,t)\) is smooth enough, so that we can differentiate this function many times as we want and each derivative is a well-defined bounded function (that has thus a fixed upper bound). Here \(L\) is the given finite length of the domain and \(T^0(x) > 0\) is a given function defining the initial condition on the temperature. The coefficient \(\kappa > 0\) is the thermal diffusivity. In general, the parameter \(\kappa\) is called the diffusion coefficient or sometimes the dispersion coefficient.
Tip
Always check the dimension of parameters that appear in the PDE. For instance, in the given heat equation, the unit of \(\kappa\) is m\(^2\)/s. (Check yourself!)
In the above case, we have imposed a Dirichlet boundary condition at \(x = 0\) and a Neumann boundary condition at \(x = L\). Since both these boundary conditions are time independent, we expect the solution to eventually reach a steady state solution \(T(x,t) \to {\tilde T}(x)\) which then remains essentially unchanged at later times. Hence, we will marching forward in time until a steady state solution is found. During this marching a transient solution will be obtained that must fulfill the initial condition.
Show that the steady state solution is \(\underset{t \to \infty}{\lim} T(x,t) = {\tilde T}(x) = 1, \quad x \in [0,1]\).
Solution
To find the stationary solution, the time derivative is set to zero. Hence, we consider the following equation
Integrating twice, we obtain the following solution
where \(a\) and \(b\) are the constants of integration. Their values can be found by means of the boundary conditions.
First, \({\tilde T}(0) = 1\), so that \(b = 1\). Second, \(\partial {\tilde T}/\partial x = 0\) at \(x = 1\), so that \(a = 0\). Hence, the steady state solution is given by
Well posedness
The above mathematical model (1.2) is well posed since we found a unique solution that is stable. An example of an ill posed problem would be the one in which the Dirichlet boundary condition at \(x = 0\) is replaced by the following homogeneous Neumann condition:
In that case, the steady state solution would be
where \(b\) remains undetermined. Hence, this solution is not unique!
Example#
Consider the BVP (1.2) and let the following parameters be given:
\(\kappa = 1\) m\(^2\)/s,
\(L = 1\) m
\(T^0 = 0\) K (uniform throughout the rod),
A function called generate_solution has been implemented that solves the heat equation numerically.
(Note that the code of this function is hidden but can be inspected by downloading this page as .ipynb-file)
The associated approach will be discussed later in this section.
Next, a plot function is defined to plot the solution to our BVP for a given time step. Run the cell below and try different time steps and observe how the spatial distribution of temperature within the rod evolves.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, fixed, widgets
# define function to generate solution to heat equation
def generate_solution():
# spatial discretization
dx = 0.1
x = np.arange(0, 1+dx, dx)
nx = len(x)
# time discretization
dt=0.005
tend=10.
nt = int(tend/dt)
# solution vector for all time steps
T = np.zeros((nt,nx))
# temperature field at t=0
T[0,:] = 0.
# loop through time
for n in range(nt-1):
# solution on the left
T[n+1,0] = 1.
# solution in the interior
for m in range(1,nx-1):
T[n+1,m] = T[n,m] + dt/dx/dx * (T[n,m+1]-2*T[n,m]+T[n,m-1])
# solution on the right
T[n+1,-1] = T[n,-1] + 2*dt/dx/dx * (T[n,-2] - T[n,-1])
return x, T
def plot_solution(x, sol, step):
fig = plt.figure()
ax = plt.axes(xlim=(0, 1), ylim=(-0.5, 1.5))
ax.plot(x, sol[min(step,1999),:], marker="o")
plt.xlabel('distance in rod [m]')
plt.ylabel('temperature [K]')
plt.title(f"temperature distribution in rod at step {step}")
plt.show()
def show_animation():
play = widgets.Play(min=0, max=400, step=1, value=0, interval=100, disabled=False)
slider = widgets.IntSlider(min=0, max=400, step=1, value=0)
widgets.jslink((play, 'value'), (slider, 'value'))
interact(plot_solution, x=fixed(x), sol=fixed(T), step=play)
return widgets.HBox([slider])
# call function to generate solution to the heat equation
# the array x defines the grid and array T contains the numerical solution
x, T = generate_solution()
# plot solution at time step 0
plot_solution(x,T,0)
# alternatively, use the slider below
# note: uncomment the following line
#interact(plot_solution, x = fixed(x), sol =fixed(T), step = widgets.IntSlider(min=0, max=400, value=0));
Domain discretization#
The first step is to discretize the domain \([0,L]\) into a finite number of grid points \(M+1\) where we intend to compute the solution of the PDE (1.2). \(M\) represents the number of spatial intervals or grid cells. We will use a uniform or equidistant grid, with a grid spacing \(\Delta x = L/M\), which is the size of each grid cell.
We will refer to one of the points in the grid as \(x_m = m\Delta x\,, m=0,\ldots,M\). The first and last grid points of the domain, \(x_0\) and \(x_M\), are called the boundary points (or in 1D the end points), while the other grid points are the internal (or inner) points. Fig. 2.17 illustrates an example of 1D uniform grid in \(x-\)space with 5 inner grid points \(x_1, \cdots, x_5\) and 2 boundary points \(x_0\) and \(x_6\). Also, there are 6 grid cells.
Likewise we discretize the time interval into a number of time steps, separated by a time increment \(\Delta t\), and compute the solution for times \(t_n = n\Delta t\,, n=1,2,3,\ldots\) until a steady state is found. Our aim is to compute the solution of the PDE for all values \(T^n_m \approx T(m\Delta x,n\Delta t)\). The function \(T^n_m\) is called the grid function and consists only of discrete values defined at grid points \(x_m\) and time steps \(t_n\). (Keep in mind that \(T(x,t)\) is continuous and smooth.)
Method of lines#
We now have a spatial grid or computational domain that approximates the physical domain \(x \in [0,L]\) and a discrete time frame. The next step it to define approximations to the partial derivatives appearing in the heat equation.
The approach to follow is to carry out the approximation of time derivatives and the approximation of spatial derivatives separately. Hence, the whole approximation process is done in two steps:
Approximate the spatial derivative first. This step is called the semi discretization step and leads to a system of first order ODEs.
The resulting system of ODEs is approximated by means of time integration. This step has been dealt with in Chapter Numerical Modelling.
This two-step process is called the method of lines (MOL)[1]. The MOL approach has a great flexibility for constructing a numerical scheme in the sense that the choice of spatial discretization and the choice of time integration are independent from each other.
A number of methods have been proposed in the field of numerical mathematics for spatial discretization. Popular methods are the finite difference, finite volume and finite element methods. In this chapter, we restrict ourselves to the traditional numerical method in hydraulic engineering, namely, the finite difference method. Another popular method is the finite element method which will be discussed in Chapter Finite Element Method.
Finite difference method#
The basic methodology of the finite difference method (FDM) is to approximate the spatial derivatives with differences of the unknowns on the grid. A variety of different approximations are possible and the choice depends on the desired accuracy (consult Chapter Numerical Modelling for further clarification).
Recall the definition of the partial derivative:
We do not have the luxury of computing the limit numerically so we must select a small but finite value for the mesh width \(\Delta x\) and approximate the derivative as
with \(T_m(t) \approx T(m\Delta x,t)\). Note that the quantity \(T_m(t)\) is continuous in time as there is no approximation in time yet! Approximation (1.3) is known as the forward finite difference scheme. This scheme is also called the one-sided approximation since \(T\) is evaluated only at \(x_{m+1} > x_m\).
We now quantify the accuracy of the forward difference scheme. By means of the Taylor series expansion we can derive the truncation error of this approximation. We expand \(T(x+\Delta x,t)\), as follows
with \(T_x = \partial T/\partial x\), \(T_{xx} = \partial^2T/\partial x^2\), etc. Here we assume that \(T(x,t)\) is sufficiently smooth. Hence, all its derivatives exist. Substitution gives
Hence, the forward difference formula is accurate only to order \(\Delta x\) and is called the first order scheme.
There are, however, two other common possibilities. The first is the backward finite difference scheme
and the second is the centred finite difference scheme or central differences
since this approximation is centred around the point of consideration \(x_m\). Central differences is an example of a two-sided approximation and we will see that, for a sufficiently small value of \(\Delta x\), this approximation leads to a more accurate numerical solution of the diffusion equation than a one-sided approximation. This does not necessarily imply that one-sided approximations are not appropriate. For instance, for advective transport, it may be appropriate to use either the forward or backward scheme. This will be discussed in Section Advection Equation.
Show that the backward and centred finite difference schemes are accurate to \(\Delta x\) and \(\Delta x^2\), respectively.
Solution
See Chapter Numerical Modelling.
The above central differences is measured using double grid size, that is, \(2\Delta x\). Alternatively, the partial derivative may also be approximated as follows
which is central differences using single grid size. However, the quantities \(T_{m\pm1/2}\) are not defined on the grid. In fact, their locations \(x_{m\pm1/2}\) are in between the grid points and are called midpoints. Hence, they must be computed by means of linear interpolation, as follows
Chapter Numerical Modelling outlined a method to find an approximation for the second derivative. Here, we will follow another approach. The approximation of the second derivative can be obtained by recalling that
Hence,
Note again the use of midpoints. Since this approximation is symmetric and centred around point \(x_m\) it is referred to as the central difference scheme for the second derivative.
Show that this central difference approximation is second order accurate.
Solution
Note that the grid function \(T_m\) represents discrete values at grid points \(x_m\), but is still continuous in time, that is, \(T_m(t)\). Furthermore, \(T_{m\pm1}(t) \approx T(x\pm\Delta x,t)\) while \(T(x,t)\) is sufficiently smooth in \(x\). Hence, we apply the Taylor series expansion to find the expansion of both \(T(x+\Delta x,t)\) and \(T(x-\Delta x,t)\). Since we are dealing with the second derivative and we also expect second order accuracy, we will expand them till the fourth order (2+2=4). Hence, we have
Plugging these expansions into the approximation yields
which proves second order accuracy.
Substituting this second order scheme into our original PDE, Eq. (1.2), we obtain
which is a semi discretization of Eq. (1.2) of the interior of the domain. Note that Eq. (1.7) is an ordinary differential equation.
Finally, we need to discretize the boundary conditions. The Dirichlet condition at first boundary point \(x_0\) is simply given by
The implementation of the Neumann condition at last boundary point \(x_M\) is less trivial. When using the central differences for spatial derivatives, we often need to consider values of the numerical solution that lie outside the computational domain in order to compute spatial derivatives on the boundaries. The usual approach is to introduce virtual points that lie outside the domain and to use both the Neumann boundary condition and the numerical scheme for the interior domain, to eliminate the values at virtual points. For example, the centred approximation for the first derivative \(\partial T/\partial x\) at boundary \(x=L\) is (substitute \(m=M\) in Eq. (1.5))
involving the virtual point \(x_{M+1} = L + \Delta x\). We also apply the semi discretization of the heat equation at the boundary point \(x_M\) (see Eq. (1.7)),
and eliminate \(T_{M+1} = T_{M-1}\) from this equation. We are left with the numerical treatment for the Neumann condition, involving only values of the numerical solution located within the domain, namely,
Summarising, we have the following semi-discrete linear system of first order ODEs with the unknowns \(T_m(t)\,, m=0,\ldots,M\),
It is standard practice to write this linear system in the matrix-vector notation
with
a vector containing \(M\) unknowns,
an \(M \times M\) discretization matrix, and
a vector with \(M\) elements. Note that vector \(\mathbf{g}\) is due to the Dirichlet boundary condition.
Symmetric matrix
You may have noticed that, apart from the last row and the last column, the discretization matrix is symmetric, that is, \(\mathbf{A} = \mathbf{A}^T\). It is possible to make this matrix completely symmetric. Recall the discretization of the homogeneous Neumann condition
This is rewritten as follows
with
If we redefine the vector with unknowns as
then the matrix of the system becomes
which is symmetric. This means that all the eigenvalues of \(\mathbf{A}\) are real. Hence, no harmonic motion shows up in the heat equation (as expected).
Using the theory of Gerschgorin we can determine the range of the eigenvalues of \(\mathbf{A}\). For a symmetric matrix \(\mathbf{A}\) with entries \(a_{ij}\) its eigenvalues are lying in the union of intervals
Thus, in our case we have
and
and we conclude that all the eigenvalues satisfy the inequality
Moreover, matrix \(\mathbf{A}\) is nonsingular, meaning that all the eigenvalues are nonzero and thus exclusively negative. As a consequence, the considered linear system of ODEs (1.8) is regarded as being damped. This property is also shared by the diffusion process because diffusion acts as a form of damping \(-\) the variation in \(c(x,t)\) reduces over time.
Recall the first exercise above. Now, we want to find the steady-state solution by means of the semi-discrete system of ODEs (1.8). The rod is 1 m long and is divided into five equally spaced cells, so that \(\Delta x = 0.2\) m. The grid consists of six points, denoted \(x_0 = 0,\cdots,x_5 = 1\). The Dirichlet boundary condition is still \({\tilde T}(x_0) = 1\) and at \(x_5=1\) the homogeneous Neumann boundary condition holds.
Now, write down the system of ODEs as given above for \(M=5\) and try to solve this system using the Numpy library, in particular the function numpy.linalg.solve.
(See also section on Boundary-value problems of Chapter Numerical Modelling.)
Alternatively, you may use the functions numpy.linalg.inv and numpy.dot.
Repeat the exercise with various values of \(\Delta x\) and also \(\kappa\). What are your conclusions?
Solution
Note that \(d\mathbf{T}/dt = 0\). Now, for the case \(M=5\) the system of equations is given by
which can be solved using a numpy.linalg function. The coding is given below.
import numpy as np
#
# n = rank of A
n=5
#
# build matrix A
A = np.zeros((n,n))
for i in range(1, n-1):
A[i, i] -= 2
A[i, i-1] += 1
A[i, i+1] += 1
A[0,0]=-2
A[0,1]=1
A[n-1,n-1]=-2
A[n-1,n-2]=2
#
# build right-hand side vector g
g = np.zeros(n)
g[0]=-1
#
# calculate inverse of A
A_inv = np.linalg.inv(A)
#
# compute solution vector T
T = np.dot(A_inv, g)
print(T)
The steady-state solution is given by
Notice that the given system of ODEs is independent of both parameters \(\Delta x\) and \(\kappa\), because of the fact that \(d\mathbf{T}/dt = 0\). Hence, the second order accurate scheme represented by the above system (1.8) for an arbitrary \(M\) provides an exact (constant) solution \({\tilde T} = b\) for the following equation
augmented with the Dirichlet boundary condition at \(x=0\) and the homogeneous Neumann condition at \(x=L\).
Time integration#
The semi discretization of the heat equation leads to a linear system of first order ODEs, (1.8). In accordance with the method of lines, the next step is to integrate this system with respect to time. Our choice for time integration would be the forward Euler scheme, since this is the most simple one. Application of forward Euler to Eq. (1.8) yields
with \(\Delta t\) the time step and \(T^n_m \approx T(m\Delta x,n\Delta t)\) a discrete grid function, which represents the wanted numerical solution to the heat equation, Eq. (1.2), at time step \(n\) and grid point \(m\). We can re-arrange this equation as follows
At each new time step, the dependent variable \(T\) at each interior grid point \(x_m\) is computed from values of \(T\) at three grid points, \(x_{m-1}\), \(x_m\) and \(x_{m+1}\), at the preceding time step. This expression defines an explicit scheme for solving the heat equation. Hence, the solution can be found at any time through a sequence of time steps. This scheme is known as the FTCS scheme, which stands for Forward in Time, Central in Space, since we have used the forward approximation in time with explicit Euler and the central differences in space.
Truncation error#
A common method to check the order of accuracy of a numerical scheme is to compute its (global) truncation error.
Definition (truncation error)
Let be given a PDE with its solution \(T(x,t)\) and let a numerical scheme for this PDE be given by
where \(L_{\Delta t, \Delta x}\) is a discrete linear operator depending on the numerical parameters \(\Delta t\) and \(\Delta x\), and \(T^n_m\) is its numerical solution at time step \(t_n = n\Delta t\) and grid point \(x_m = m\Delta x\).
The truncation error of the scheme is defined as
Let us consider the FTCS scheme (1.11). We can derive its truncation error, as follows
with \(T(x,t)\) a sufficiently smooth, exact solution to the heat equation (1.2). Because of this, we can use the Taylor series expansions of \(T(x_m,t_{n+1})\) and \(T(x_{m\pm1},t_n)\) to further evaluate the truncation error. So,
and, thus
Hence, the FTCS scheme is first order in \(\Delta t\) and second order in \(\Delta x\). In practice, this means we can discretize the spatial domain quite coarsely but we must discretize the time interval more finely in order to achieve a required accuracy.
We consider the ODE related to the Neumann boundary condition at \(x=L\), which is given by (see (1.8))
This ODE can be approximated using the forward Euler scheme which yields
Compute the associated truncation error.
Solution
First, we replace the numerical solution \(T^n_m\) by the exact solution \(T(x_m,t_n)\). Hence, we have
The exact solution is supposed to be smooth so that we can use the Taylor series expansions for the functions \(T(L,t+\Delta t)\) and \(T(L-\Delta x,t)\). We expect a first order accuracy in time and hence the Taylor series of \(T(L,t+\Delta t)\) is expanded till the second order term (”first derivative + first order accuracy”). With respect to the Taylor series of \(T(L-\Delta x,t)\), we will expand this till the third order (=1+2) because we expect at most second order accuracy (it can be less!). Hence,
and
Substitution gives
where we have used the homogeneous Neumann condition at \(x=L\), namely, \(T_x\,(L,t) = 0\).
So, we conclude that the numerical scheme related to the Neumann condition at \(x=L\) is first order accurate both in time (as expected!) and in space (not expected!).
Our first thought might be to conclude that our overall approximation in the whole(!) domain, including the boundary points, might have lost second order accuracy. However, this is not the case. Generally, we achieve second order accuracy in the overall numerical method even if the truncation error is \(\mathcal{O}(\Delta x)\) at a single boundary point as long as it is \(\mathcal{O}(\Delta x^2)\) everywhere else. So, we can often get away with one order of accuracy lower in the local error for the boundary conditions than what we have elsewhere.
Consistency and convergence#
Definition (consistency)
A numerical scheme is called consistent if and only if
According to this definition, the FTCS scheme (1.11) is consistent with the heat equation (1.2).
The above error analysis shows how well the FTCS scheme approximates the original heat equation in the limit of \(\Delta x, \Delta t \to 0\). However, it does not show how well the numerical solution approximates the exact solution of the heat equation. This is a matter of convergence.
Definition (convergence)
A numerical scheme as described above is said to be convergent if
In the study of convergence of a numerical scheme, there are two issues to be considered.
Firstly, whether the numerical approximation is consistent with the PDE we wish to solve.
Secondly, whether the considered numerical scheme is stable, that is, its numerical solution remains bounded in case of endless repeating of this scheme, with fixed \(\Delta t\) and \(\Delta x\). This topic of stability will be investigated in some detail in the next section.
With these two considerations, we can apply the following
Basic rule
Since a proof of convergence is the hard one, effort in the other two aspects is easily elaborated. Consistency defines a relation between the numerical scheme and the PDE we wish to solve, whereas stability establishes a relation between the computed solution and the exact solution of the numerical scheme. It should be emphasised that stability is a requirement just on the numerical scheme and does not involve any condition on the PDE.
Stability#
We may write the FTCS scheme (1.11) and (1.13) in the matrix-vector notation, as follows (see also (1.9))
with \(\mathbf{T}^n = (T^n_1,\cdots,T^n_M)^{\text{T}}\) being the vector representing the numerical solution at time step \(n\) and \(\mathbf{A}\) is the discretization matrix as given by Eq. (1.10). Matrix \(\mathbf{I}\) is the identity matrix. A way to check stability is to consider the eigenvalues \(\lambda_m\) of the matrix \(\mathbf{I} + \Delta t \mathbf{A}\) and subsequently require \(|\lambda_m| \leq 1\,, \forall m\). However, the matrices \(\mathbf{A}\) and \(\mathbf{I}\) are an \(M \times M\) matrix with \(M = L/\Delta x\) so that its dimension is growing as \(\Delta x \to 0\). Hence, we might not be able to determine eigenvalues of such a large matrix, because of the cumbersome and laborious algebra.
Derive stability condition using eigenvalues
For our discretization matrix \(A\) it is still possible to derive a stability condition. Specifically, using the Gerschgorin’s theorem, we found the range of eigenvalues of \(A\) which is given by
The smallest eigenvalue equals \(-4\kappa/\Delta x^2\). Hence, the system of equations is stable under the condition
The last inequality is always true while the first inequality is true if
From a physical point of view, the solution to the heat equation (1.2) cannot be negative, that is, we must have \(T(x,t) \geq 0\, , \forall x \in [0,L]\,, \forall t \geq 0\). This implies that our finite difference scheme should share the same property, so that the numerical solution remains non-negative as time progresses. As a consequence, spurious oscillations (wiggles) will not occur. It is well known that in the case of diffusion processes, wiggles in the solution will blow up and thus will be unbounded in finite time. To prevent this instability, we have to require that our numerical scheme should not exhibit spurious oscillations.
Recall (1.12). We can rewrite this equation as follows
Assume that \(T^n_m \geq 0\) for all grid points \(m = 1,\ldots,M-1\), then the last two terms on the right-hand side are non-negative. Furthermore, if
then \(T^{n+1}_m \geq 0\) for all grid points \(m = 1,\ldots,M-1\). Hence, by induction it follows that the numerical solution is non-negative at all times.
This resulted stability condition is
and states that for a given \(\Delta x\), the allowed value of \(\Delta t\) must be small enough to keep the FTCS scheme (1.11) stable. For this reason, the FTCS scheme is said to be conditionally stable.
Verify the above stability condition for grid point \(m=M\).
Solution
We rewrite (1.13) as follows
or
We assume that the values \(T^n_M\) and \(T^n_{M-1}\) are non-negative. While the second term on the right-hand side is non-negative, the first term is non-negative as long as the following condition holds
which is the same condition for the inner points \(m = 1,\ldots,M-1\).
In this exercise we check the stability of the FTCS scheme for our BVP (1.2) with the following parameters:
\(\kappa = 1\) m\(^2\)/s
\(L = 1\) m
Furthermore, the rod is divided into 10 grid cells (\(M=10\)).
What is the maximum time step \(\Delta t\) to get a stable solution?
Feel free to use the widget below by adjusting the time step interactively. (Note: the plot below shows the temperature distribution in the rod at last time step which is the steady state solution \({\tilde T}(x) = 1\) K.)
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, widgets
# define function to carry out the FTCS scheme
def ftcs(dt):
# spatial discretization
dx = 0.1
x = np.arange(0, 1+dx, dx)
nx = len(x)
# time discretization
tend=3.
nt = int(tend/dt)
# solution vector for all time steps
T = np.zeros((nt,nx))
Texact = np.zeros(nx)+1.
# temperature field at t=0
T[0,:] = 0.
# loop through time
for n in range(nt-1):
# solution on the left
T[n+1,0] = 1.
# solution in the interior
for m in range(1,nx-1):
T[n+1,m] = T[n,m] + dt/dx/dx * (T[n,m+1]-2*T[n,m]+T[n,m-1])
# solution on the right
T[n+1,-1] = T[n,-1] + 2*dt/dx/dx * (T[n,-2] - T[n,-1])
fig = plt.figure()
ax = plt.axes(xlim=(0, 1),ylim=(np.ceil(T[-1,:].min()) - 0.1, np.ceil(T[-1,:].max()) + 0.1))
ax.plot(x,T[-1,:],label='numerical solution', color="red", marker="o")
ax.plot(x,Texact, label="exact solution", color="black")
ax.legend()
plt.xlabel('distance in rod [m]')
plt.ylabel('temperature [K] at last step')
plt.title(f"FTCS scheme with dt = {dt} s")
plt.show()
# create an interactive slider
interact(
ftcs,
dt=widgets.FloatSlider(value=0.003, min=0.001, max=0.006, step=0.0001, description="dt",readout_format='.4f'),
);
Solution
The stability condition of the FTCS scheme reads
Furtheremore, \(M = 10\) and \(L = 1\) m, and so \(\Delta x = 0.1\) m. With \(\kappa = 1\) m\(^2\)/s, we have