Week 1: Finite differences and Euler’s method
\[ \renewcommand{\vec}[1]{\boldsymbol{#1}} \newcommand{\td}[2]{\frac{\mathrm{d}#1}{\mathrm{d}#2}} \newcommand{\tdd}[2]{\frac{\mathrm{d}^2#1}{\mathrm{d}#2^2}} \newcommand{\pd}[2]{\frac{\partial#1}{\partial#2}} \newcommand{\pdd}[2]{\frac{\partial^2#1}{\partial#2^2}} \]
Overview
To start off the unit, we dive into finite differences. This is a general method for approximating derivatives of mathematical functions with algebraic formulae. This will be a crucial step when it comes to deriving numerical methods for solving ODEs and PDEs. We also introduce the idea of truncation error, explain why this is useful, and how it can be calculated using Taylor expansions. Finally, we recap Euler’s method, a simple yet useful approach for numerically approximating the solution of time-dependent ODEs.
Supplementary material
Exercises
Consider the following finite-difference approximation of the second derivative \[ \left.\tdd{u}{x}\right|_{x=a} = \frac{2 u(a) + C u(a+\Delta x) + 4u(a + 2 \Delta x) - u(a + 3\Delta x)}{(\Delta x)^2} \] First, take \(C = 5\). By calculating the truncation error analytically or numerically, determine whether this is a valid approximation formula. If the formula is valid, determine the order of the approximation. Is the approximation formula valid when \(C = -5\)? If so, what is the order of the approximation?
Use Euler’s method to solve the ODE system given by \[ \td{u_1}{t} = -\omega u_2, \qquad \td{u_2}{t} = \omega u_1, \tag{1} \] over the time range \(0 \leq t \leq 2 \pi\). The initial condition can be set to \(u_1(0) = 1\) and \(u_2(0) = 0\). The exact solution to this problem is given by \(u_1(t) = \cos(\omega t)\) and \(u_2(t) = \sin (\omega t)\). You can use this exact solution to validate your code and check that it is working correctly.
- Set \(\omega = 1\) and use \(N_t = 100\) time steps to find a numerical solution. Plot \(u_1\) as a function of time using your numerical solution and the exact solution.
- Keeping \(\omega = 1\), create similar plots when \(N_t\) is decreased to 50 and increased to 1000. From the plots, what happens to the accuracy of the numerical solution as \(N_t\) increases?
- Set \(\omega = 1\) and use \(N_t = 100\) time steps to find a numerical solution. Plot \(u_1\) as a function of time using your numerical solution and the exact solution.
For the ODE system given by (1), use numerical experimentation to determine how different values of \(\omega\) impact the solution. If \(u_1\) represents the displacement of a mass on a spring, then what does \(\omega\) physically correspond to?
For the ODE system given by (1), compute the error between your numerical solution and the exact solution at time \(t = 2 \pi\); that is, calculate \(\varepsilon = |u_1^{N_t} - 1|\), where \(u_1^{N_t}\) is the numerical solution at time step \(N_t\).
- For a fixed value of \(\omega = 1\), how fast does the error decrease with \(\Delta t\)? Is it linear, quadratic, etc? Is this what you expect based on the truncation error of Euler’s method?
- For a fixed number of time steps, \(N_t = 10^4\), how does the error change as \(\omega\) increases? Can you provide a justification of the results you see? For this problem, try setting \(\omega = 1, 2, 4, 8, 16\).
- For a fixed value of \(\omega = 1\), how fast does the error decrease with \(\Delta t\)? Is it linear, quadratic, etc? Is this what you expect based on the truncation error of Euler’s method?
Coursework-style problem
Consider a pack of \(N\) identical lithium-ion batteries that are connected in series. Let \(i = 1\) denote the first battery in the pack and \(i = N\) denote the last battery in the pack. The other batteries are labelled with \(i = 2, 3, \ldots, N - 1\). Due to the geometry of the pack, only batteries 1 and \(N\) are in contact with the environment. Battery \(i\) is in contact with batteries \(i-1\) and \(i + 1\).
Each battery in the pack acts as a resistor. When an electrical current is drawn from the battery pack, each battery will heat up due to Joule heating. The heat that is generated in battery \(i\) will be transferred to batteries \(i-1\) and \(i + 1\), causing these two batteries to heat up even more. Heat is transferred from one battery to the next until it reaches batteries 1 and \(N\), which can then transfer the heat to the air. The more batteries there are, the longer the generated heat remains in the battery pack, and the hotter the pack will get. If the temperature of any battery exceeds 60 \(^\circ\)C, then it will start to degrade. Battery degradation is extremely dangerous, as this can lead to the release of toxic and explosive gases.
The aim of this exercise is to simulate the temperature of a battery pack in order to determine whether it is expected to be safe to use.
Assume that a constant current \(I\) is being drawn from the battery pack. If \(T_i(t)\) denotes the temperature of battery \(i\) at time \(t\), then the temperature of each of the batteries can be obtained by solving the ODE system given by \[\begin{align} c \td{T_1}{t} &= I^2 R - h(T_{1} - T_\text{air}) - h(T_{1} - T_{2}), \tag{2} \\ c \td{T_i}{t} &= I^2 R - h(T_{i} - T_{i-1}) - h(T_{i} - T_{i+1}), \quad i = 2, 3, \ldots, N - 1, \tag{3} \\ c \td{T_N}{t} &= I^2 R - h(T_{N} - T_\text{air}) - h(T_{N} - T_{N-1}). \tag{4} \end{align}\] Here, \(c\) is the specific heat of the battery, \(R\) is the resistance of the battery, \(h\) is the heat transfer coefficient, and \(T_\text{air}\) is the temperature of the air. Values for these parameters are given in the table below. You can assume that the initial temperature of each battery is given by \(T_\text{air}\); that is, \(T_i(0) = T_\text{air}\).
| \(I\) (A) | \(c\) (J/\(^\circ\)C) | \(h\) (W/\(^\circ\)C) | \(T_\text{air}\) (\(^\circ\)C) | \(R\) (Ohms) |
|---|---|---|---|---|
| 3 | 300 | 0.7 | 23 | 0.4 |
In the exercises below, assume that the current is being drawn for one hour.
Consider a battery pack with \(N = 5\) batteries.
- Numerically solve the ODE system in (2)–(4). Create a single plot that shows the temperature of each battery as a function of time.
- At each time \(t\), compute the maximum temperature across all 6 batteries. Plot the maximum temperature as a function of time.
- Use your results to determine whether the battery is safe to use. Answer: yes.
Is a battery pack with \(N = 50\) batteries safe to use? Answer: no.
What is the largest number of batteries that a pack can contain yet still be safe to use? Answer: 6.
Answers to selected exercises
\(C = 5\) does not lead to a valid differentiation formula. \(C = -5\) does lead to a valid differentiation formula that is second-order accurate, with \(TE = O((\Delta x)^2)\).
- Use the exact solution to validate your code and see what the numerical solution should look like.
- The accuracy of the solution increases as \(N_t\) increases.
The parameter \(\omega\) corresponds to the natural frequency of the mass-spring system.
- The error decrease linearly with \(\Delta t\).
- The error increases quadratically with \(\omega\). This result is agrees with theory, as the truncation error is proportional to \(\ddot{u}_1\). From the exact solution, \(\ddot{u}_1\) is proportional to \(\omega^2\). This result implies that the accuracy of Euler’s method not only depends on the size of \(\Delta t\), but also how fast the derivative of the solution is changing.
- The error decrease linearly with \(\Delta t\).