Ordinary differential equations

I. ODEs as Mathematical Models

Mathematical modeling is the process of developing mathematical models (equations) that describe physical phenomena or other systems. A mathematical model help to study, understand, and explain a complex system and to make predictions about behaviours of the system (both in the past and future). To a large extent the success of mathematical modeling is due to the fields of numerical analysis, scientific computing, and more recently uncertainty quantification. These fields aim at: 1) obtaining valid mathematical models; 2) translating the models into accurate and efficient algorithms (referred to as computational models); and 3) implementing the algorithms on computers. The computer outputs will then enable scientists to produce predictions about real-world processes and phenomena. The reliablility and feasibility of the predictions depend on: i) the validity of mathematical models; ii) the ingenuity of computational models; and iii) the power of the computer used to produce the predictions. It is important to note that the feasibility of predictions depends on both the power of computers and the ingenuity of the numerical methods used in the algorithms. For example, solving an \(n \times n\) linear system of equations requires \(\sim n^3\) floating point operations (FLOP). If we use a modern computer with a performance on the order of \(10^{12}\) floating point operations per second (FLOPS), then solving a system with \(n=10^6\) unknowns would take \(10^{18} \, / \, 10^{12} \, / \, 60 \, / \, 60 \, / \, 24 \, / \, 365 \, \approx \, 1\) year. On the other hand, if we use the incredibly ingenious multigrid method which only requires about \(C\times n\) operations to solve the same system it would instead take about \(C\) millionth of a second! This kind of reduction is crucial in reducing the use of computing resources and the time to obtain a prediction. The latter has led to much shorter design and inovation cycles in engineering and production.

The process of obtaining a valid mathematical model for a given system starts with applying a set of rules or laws that are thought to be valid for the system under study. For example when modeling a fluid flow the conservation of mass is such a rule, or when modeling the vibration of a guitar string we apply the second law of motion. In some cases, for instance in economics, the rules may not be motivated by physics. A simple rule that you use every day is the conservation of money. Suppose you have $5 to spend, and you go to store to buy apple, which costs $2 per lbs. How much apple can you buy? You formulate the problem: x lbs of apple costs 2x and I have $5. Then conservation of money states that 2x=5. This is a mathematical model (an equation). Solving your model, you realize that you can buy 5/1.89 = 2.5 lbs of apple.

For more complicated problems the models are more complicated and it may not always be clear what rules/laws one should choose to enforce. Then the modeling process could become iterative:

1. Choose an admissible mathematical model.
2. Solve the model by an accurate and efficient computational model on a computer.
3. Compare the computer output to experimental data.
4. If the simulated result is “close” to the experimental data, then the selected model is valid. Otherwise, go back to step 1 and try to improve the model.

In practice the above steps need to be done in a systematic way through various processes, including calibration, validation, verification, and prediction.

Differential equations (both ordinary and partial) are important classes of mathematical models. In this section, we focus on ordinary differential equations (ODEs). ODEs serve as mathematical models to many systems. Examples include prey-predator systems, carbon dating, chemical reactions, etc. In Homework5 you will experiment with some different models for the group dynamics of a flock of birds. In addition to serving as mathematical models, ODEs also often appear when solving partial differential equations (PDEs). In reality it is not possible to solve (complicated) ODEs by analytical methods. In such cases we will need to find approximate solutions to ODEs by numerical methods.

II. Numerical methods for solving ODEs

Consider the following initial value problem, consisting of a first-order ODE and an initial condition:

\begin{eqnarray} \left\{ \begin{array}{l l} {\text{ODE:}} & \, \, u' = f(t,u(t)), \hskip 1.4cm t \in [0,T]\\ {\text{IC:}} & \, \, u(0) = u_0. \end{array} \right. \end{eqnarray}

We first discretize the interval \([0,T]\) into \(n+1\) equally spaced grid points \(\{ t_i = i \, \Delta t \}_{i=0}^n\), where \(\Delta t = T/n\) is the time step. We then use Taylor’s formula and expand \(u(t+\Delta t)\) around \(t\), using only two terms

$$u(t+\Delta t) = u(t) + \Delta t \, u'(t) + \mathcal{O} (\Delta t^2).$$

When \(\Delta t\) is small, or equivalently \(n\) is large, the term \(\mathcal{O} (\Delta t^2)\) is negligible, and hence we can use this formula to approximate the time derivative on the left hand side of the ODE:

$$u'(t) \approx \frac{u(t+\Delta t)-u(t)}{\Delta t}.$$

Now let \(u_i\) denote the approximation of \(u(t_i)\), that is, let \(u_i \approx u(t_i)\). Given \(u_0 = u(t_0)\) by the initial condition, we can approximate \(u(t)\) at all grid points \(\{t_1, t_2, \dotsc, t_n \}\) on the interval \([0,T]\) by:

\begin{eqnarray*} u_0 &=& \text{given by IC}, \\ u_{i+1} &=& u_i + \Delta t \, f(t_i,u_i), \ \ i= 0,1,\ldots, n-1 \end{eqnarray*}

The above method is called forward Euler. It is the simplest method for solving ODEs.

There are a veriety of numerical methods for solving ODEs, such as Runge-Kutta methods and multistep methods. You can get started with reading my lecture notes to learn more about Runge-Kutta methods. For further reading I would recommend this and this text books.

Accuracy

An important aspect of a numerical method is its accuracy. We want the error in the approximation computed by a numerical method to be small. In the case of forward Euler, it turns out that the error is proportional to \(\Delta t\). In other words the error decreases linearly with respect to \(\Delta t\). We say that the “order of accuracy” of forward Euler is one.

This might be surprising, because Taylor’s formula tells us that the error in the approximation over one time step is proportional to \((\Delta t)^2\). What relly happens is that we compute the solution tp to some final time \(T\). This means we need to take \(n = \frac{T}{\Delta t}\) time steps. If we commit an error of order \((\Delta t)^2\) in each time step, the final error at the final time \(T\) will be proportional to \(n \, (\Delta t)^2 \sim T \, \Delta t\). Hence, the overal order of accuracy is one. This is called the accumulation of error.

In general, a numerical method for solving ODEs is said to be \(p\)-th order accurate if the approximation error is proportional to \((\Delta t)^p\), with \(p \ge 1\). The larger the \(p\), the more accurate the moethod. Assuming there is a maximum allowable error (denoted by \(\varepsilon_{\text{TOL}}\)), we obtain an upper limit on \(\Delta t\). In other words, to achieve the given tolerance, the time step cannot be larger than this upper limit. This is called the accuracy limit.

Stability

Stability is another important feature of a numerical method. We want a numerical method to be stable. Stability is a property that together with accuracy will gurantee convergence. Convergence would ensure that the numerical solution will converge to the correct solution as \(\Delta t \rightarrow 0\). Without stability, the approximation error may drastically grow, causing an exponential growth in the deviation of the approximate solution from the exact solution. We refer to this as the explosion of the numerical error. Stability may impose a stability limit on \(\Delta t\) (in addition to the accuracy limit). To illustrate the stability limit we consider the following two exercises.

Exercise 1. Consider the Dahlquist model problem \(u'(t) = - \lambda u(t), \ \ \lambda > 0\) and \(u_0 = 1\).

1. Find the analytical solution to the problem.
2. What is the steady state? Is the solution of the continuous problem stable?
3. Use Euler’s method to compute the approximate solution at the first few time steps, for instance for \(\lambda = 30\) and with \(\Delta t = 0.1\). Do you observe any strange behavior?
4. Derive the stability limit on \(\Delta t\) by enforcing \(G_i = | \frac{u_{i+1}}{u_i}| < 1\). Can you figure out what this condition mean?

Try to do this exercise yourself, and then look at my lecture notes to compare your results with mine.

Exercise 2. Consider another numerical method, which is pretty similar to the forward Euler:

\begin{eqnarray*} u_0 &=& \text{given by IC}, \\ u_{i+1} &=& u_i + \Delta t \, f(t_{i+1}, u_{i+1}), \ \ i= 0,1,\ldots, n-1 \end{eqnarray*}

This method is called backward Euler.

1. Try to derive the above formula using Taylor’s expansion.
2. Repeat the stability analysis (steps 1–3 in the above exercise) for backward Euler on the Dahlquist model problem, i.e. with \(f(t,u(t)) = - \lambda \, u(t)\).
3. What is the major difference between forward Euler and backward Euler in terms of stability? When do you think this method can be useful? What is the additional complication compared to forward Euler when applied to an ODE with a more complicated right hand side?