Numerical Analysis Project in ODEs for Undergraduate Students

Abstract

Designing good projects involving programming in numerical analysis for large groups of students with different backgrounds is a challenging task. The assignment has to be manageable for the average student, but to additionally inspire the better students it is preferable that it has some depth and leads to them to think about the subject. We describe a project that was assigned to the students of an introductory Numerical Analysis course at the University of Iceland. The assignment is to numerically compute the length of solution trajectories of a system of ordinary differential equations with a stable equilibrium point. While not difficult to do, the results are somewhat surprising and got the better students to get interested in what was happening. We describe the project, its solution using Matlab, and the underlying mathematics in some detail. Further, we discuss the pedagogical aspects of the project and the results in terms of its success and shortcomings.

1 Background

The project we describe in this paper was assigned to the undergraduate students of Numerical Analysis, a 6 ECTS unit course at the University of Iceland with approximately 150 students. The responsibility for this course is within the Faculty of Physical Sciences in the School of Engineering and Natural Sciences and it is mandatory for all BSc. students of Mechanical-, Industrial-, Chemical-, and Civil and Environmental Engineering as well as for all students of Physics, Engineering Physics, Mathematics, Applied Mathematics, and Mathematics and Mathematical Education. It is an elective course for students of Geophysics, Computer Science, Software Engineering, Electrical and Computer Engineering, and Chemistry.

As can be seen from this long lists the preparation and interests of the enrolled students vary considerably. Usually, the students enroll in the course in the fourth semester of six in total to complete a bachelors degree. Prerequisites for Numerical Analysis are one course in Computer Science (programming in Matlab or Python), one course in Linear Algebra, and two courses (three recommended) in Calculus. All of the students are thus familiar with applying linear algebra and calculus, but the mathematics students have also studied the theoretical aspects of these disciplines in some detail in the framework of metric spaces, ring theory, etc. Students of Mathematics and Applied Mathematics are required to enroll simultaneously in the course Theoretical Numerical Analysis (2 ECTS units), where rigid mathematical proofs of most of the material covered in the Numerical Analysis course are studied.

The Course Description is:

Fundamental concepts on approximation and error estimates. Solutions of systems of linear and non-linear equations. PLU decomposition. Interpolating polynomials, spline interpolation and regression. Numerical differentiation and integration. Extrapolation. Numerical solutions of initial value problems of systems of ordinary differential equations. Multistep methods. Numerical solutions to boundary value problems for ordinary differential equations.

and the Learning Outcomes are:

Knowledge and understanding: To complete this course the student should be able to

• 1. define, explain and give examples of the main concepts of the course, such as error, matrix factorization, interpolating polynomial, and finite differences,

• 2. state and explain the main results of the course, for example by stating Newton’s Method and the Secant Method and estimate the errors, state algorithms to solve boundary value problems using finite differences and verify the degree of the approximation.

Skills: To complete this course the student should be able to

• 1. formulate a simple mathematical problem as a numerical problem, implement it on a computer, and compute an approximate solution,

• 2. estimate the error of numerical solutions,

• 3. use computer software, such as the Anaconda Python platform or Matlab, for programming, computing, and performing numerical experiments,

• 4. validate the results of numerical computations,

• 5. use the concepts and the results of the course to develop and advance algorithms for simple problems the student has not seen before.

In the course two larger group projects count for 30% of the final grade. The assignment we describe here was the third and last part of the second project. The other two parts were to

• 1. implement adaptive integration using the trapezoidal- and the Simpson’s rule and test it for some integrals and

• 2. implement the shooting method and use it to solve a few boundary value problems.

Adaptive integration and the shooting method are discussed in sufficient detail in the lectures and in the textbook used in the class [20] to make them rather easy to do.

2 The Project and Its Solution

In the project the time-reversed van der Pol oscillator

(1)

is considered. It has an exponentially stable equilibrium at the origin and an unstable periodic orbit around it; see below for more details. The project has three objectives, which we describe below together with its solution using Matlab.

2.1 Objective I

The first objective of the assignment was to analyze the system (1) by drawing solution trajectories, both forward and backwards in time. This is easily achieved by first defining in f.m

and then typing in the command window, here using the initial value $(2,0)^T$ for $\mathbf{x}'=\mathbf{f}(\mathbf{x})$ and $(1,1)^T$ for $\mathbf{x}'=-\mathbf{f}(\mathbf{x})$ and in both cases integrating over the time-interval [0, 20].

The plots produced are drawn in Fig. 1. A few comments are in order. Usually, one defines the function f as a function of both time t and space x, even though the system is autonomous, i.e. $\mathbf{f}$ does not depend explicitly on time. We follow this tradition here, for otherwise we could not use the Matlab solver ode45 directly. The system $\mathbf{x}'=-\mathbf{f}(\mathbf{x})$ is a time-reversion of the time-reversed van der Pol oscillator from (1), i.e. it is the van der Pol oscillator. It is well known that the van der Pol oscillator has a stable periodic orbit and it can be clearly seen in Fig. 1 (right). Since the periodic orbit is stable all “normal” initial values, except the unstable equilibrium point $(0,0)^T$, will converge to the orbit. The system $\mathbf{x}'=\mathbf{f}(\mathbf{x})$ is the time-reversed van der Pol oscillator and the stable periodic orbit of the van der Pol oscillator becomes an unstable orbit and the unstable equilibrium point at the origin becomes stable. Thus, any initial value inside of the periodic orbit will converge to the equilibrium and outside of the orbit will diverge. Note that the Matlab solver ode45 (and others) reports problems, i.e. “Unable to meet integration tolerances without reducing the step size below the smallest value allowed” in the integration for various initial values and time-intervals. This can be overcome by using a solver with fixed time-steps, to be delivered in the next objective.

Fig. 1.

Left: A solution trajectory of system $\mathbf{x}'=\mathbf{f}(\mathbf{x})$ with $\mathbf{f}$ from (1), starting at $(2,0)^T$ and integrated numerically over the time-interval [0, 20] using Matlab’s ode45. Right: A solution trajectory of system $\mathbf{x}'=-\mathbf{f}(\mathbf{x})$ with $\mathbf{f}$ from (1), starting at $(1,1)^T$ and integrated numerically over the time-interval [0, 20] using Matlab’s ode45.

2.2 Objective II

The second objective of the project was to program the standard Runge-Kutta method of order 4, commonly abbreviated RK4. In more detail, an implementation for the function

should be given, where f is the function $\mathbf{f}$ in $\mathbf{x}'=\mathbf{f}(t,\mathbf{x})$, xi is the initial value at time a, and [a,b] is the time-interval over which the solution is computed, and n is the number of time-steps used. The output [t,w] should be the transposes of the outputs of the Matlab solver ode45 (column vector notation). This is an easy task whose solution is given in the Appendix in Program I. Using RK4 to plot solution trajectories as in Objective I now becomes, using 500 time-steps:

This delivers trajectories comparable to the ones in Fig. 1. Note, that much fewer time-steps, e.g. 100, results in much less accurate results, see Fig. 2. We now move to the main objective of the problem.

Fig. 2.

Left: A solution trajectory of system $\mathbf{x}'=\mathbf{f}(\mathbf{x})$ with $\mathbf{f}$ from (1), starting at $(2,0)^T$ and integrated numerically over the time-interval [0, 20] using RK4 and 100 time-steps (step-size 20/100 = 0.2). Right: A solution trajectory of system $\mathbf{x}'=-\mathbf{f}(\mathbf{x})$ with $\mathbf{f}$ from (1), starting at $(1,1)^T$ and integrated numerically over the time-interval [0, 20] again using RK4 and 100 time-steps (step-size 20/100 = 0.2).

2.3 Objective III

On a uniform $101\times 101$ grid on $[-3,3]\times [-8,8]$, compute the length of the solution trajectories to (1) integrated over a time-interval of length 4 and using RK4 with 100 time-steps. Note that 100 time-steps for a time-interval of length 4 corresponds to the 500 time-steps for a time-interval of length 20 used above, because both have step-size 0.04. The solution to system (1), starting at at time zero, is a function

(2)

for all t in the definition domain of (dependant of either $\mathbb {R}$ or the maximum domain before becomes infinite). The length of the trajectory on the time-interval [0, T] is well known to be defined as

where $\Vert \cdot \Vert$ denotes the Euclidian norm on $\mathbb {R}^2$. Because of (2) we can substitute for and the formula for becomes

(3)

Since [t,w]=RK4(@f,xi,0,T,n) delivers in its ith column, where and $t_i=(i-1)$T/n for $i=1,2,\ldots ,$n + 1, the integral in (3) can easily be approximated using numerical integration. In the project it was suggested to use the composite Simpson’s Rule

$$\begin{aligned} \int _0^T g(\tau )d\tau \approx \frac{h}{3}\left( g(t_1)+g(t_{n+1})+4\sum _{i=1}^{n/2}g(t_{2i})+2\sum _{i=1}^{n/2-1}g(t_{2i+1})\right) , \end{aligned}$$

where $h=T/n$ is the length of the time-steps. Note that for this formula we need that n is an even number. Since many of the trajectories will be very long, indeed so long that the numerical solver will fail to deliver a numerical value, it is useful to substitute e.g.  (Matlab’s min interprets NaN as larger than any number). The implementation is now simple and is given in the Appendix in Program II. Typing

Fig. 3.

The length of the trajectories of system (1) plotted as a function of starting position (x, y).

in the command window, where TraLengths is the program from the Appendix, f is the function from f.m, $[-3,3]\times [-8,8]$ is the the area in the plane where a function is computed and plotted, T = 4 is the interval of integration in (3), and n=100 is the number of time-steps used in the numerical integration, now delivers after a short time (12.5 s on a computer with an i7-7700K CPU) Fig. 3. After having produced this figure the less interested students have finished the project. The more interested ones, however, now might start to ask themselves and the instructor why this function looks like it does? Indeed, many of them did ask the instructor.

3 Study of the Results

Let us discuss the results from Objective III and the function V computed in more detail.

Because the origin is a stable equilibrium, solutions slow down close to it. Indeed, they become so slow that the limit when $t\rightarrow \infty$ is never reached. The reason for this is that since $\mathbf{f}$ is continuous and , $\mathbf{f}$ and then also $\mathbf{x}'$ are small close to the origin. Therefore trajectories that start close to the origin are short and .

The function V is a so-called Lyapunov function for the system. This means that it has a local minimum at the stable equilibrium at the origin and that it is decreasing along all solution trajectories in a neighbourhood of the equilibrium. The local minimum is intuitively clear and the decrease can be seen from the following calculations:

(4)

and because solution trajectories are unique, i.e. , we get

Thus, by (4) and the Fundamental Theorem of Calculus,

Since for all in a neighbourhood of the origin, we have for all such that for sufficiently large $T>0$ and . Thus

and V is decreasing along solution trajectories in a neighbourhood of the stable equilibrium at the origin, i.e. the mapping is decreasing.

The stability theory of Lyapunov and Lyapunov functions are covered in most textbooks on dynamical systems and/or control theory, e.g. [16, 19, 21], and the interested student can be pointed to them for additional information. For an attractor, like our stable equilibrium at the origin for the time-reversed van der Pol system, the sublevel sets of a Lyapunov function serve as “traps”. Once inside the component of a sublevel set, that includes the origin and does not extend to the boundary of the domain of the Lyapunov function, the solution cannot escape. This comes because the solution is decreasing along solution trajectories and therefore cannot climb over the edge/boundary of the sublevel set.

The theory of complete Lyapunov functions even tells us that every system given by an ODE possesses a complete Lyapunov function that goes a long way in characterizing the qualitative behaviour of the system. Indeed, this holds true for very general dynamical systems. A complete Lyapunov function is a scalar-values function from the whole state-space that is non-increasing along along all solution trajectories and strictly decreasing where possible. Note that, e.g. for a periodic orbit, it cannot be strictly decreasing. In general it is strictly decreasing along all solution trajectories on the part of the state-space where the flow is gradient-like and constant on every transitive component of the chain-recurrent set. The theory of complete Lyapunov functions was developed by Auslander, Conley, and Hurley [1, 6, 13,14,15], see also [18].

Computing Lyapunov functions by using results from converse theorems in dynamical systems, i.e. theorems guarantying the existence of certain kinds of Lyapunov functions for systems with particular stability properties, using integrals or sums over solutions trajectories has been studied in numerous publications [2,3,4,5, 7,8,9,10,11,12, 17]. A central issue is the verification of the properties of a Lyapunov function for the function computed. One way to do this is to interpolate the values computed over the simplices of a triangulation. Essentially, one demands $\nabla V(\mathbf{x})\cdot \mathbf{f}(\mathbf{x}) + E_\mathbf{x}\Vert \nabla V(\mathbf{x})\Vert _1 < 0$ at all vertices of the triangulation, where the error $E_\mathbf{x}\ge 0$ assures that not only

at the vertices, but for any $\mathbf{x}$ in the domain of the triangulation. Here $\Vert \mathbf{x}\Vert _1=|x_1|+|x_2|$ and the function V is defined on the simplices by using convex interpolation of the values at the vertices over the whole simplex. A detailed discussion of this method is beyond the scope of this paper, but the interested reader can have a look at [11, 12].

By adding to the code for the function TraLengths as described in Program III in the Appendix, a verification of the interpolation of the values as described in [11, 12] is carried out. If one want to implement this for a different two-dimensional system than (1), only the function f in f.m and the function B in Program III have to be modified accordingly. The modification of f is obvious because it is the right-hand side of (1), and the return value for B(r,s,xm,ym) should be an upper bound, preferably close, on

$$\begin{aligned} \max _{i=1,2 \atop |x_1|\le \mathtt{xm},|x_2|\le \mathtt{ym}}\left| \frac{\partial ^2 f_i}{\partial x_\mathtt{r}\partial x_\mathtt{s}}(x_1,x_2)\right| . \end{aligned}$$

The results of this verification for our system (1) with the same parameters as in Objective III is plotted in Fig. 4 (left) and with a higher spatial resolution, i.e. xres=301; yres=301;, on the right.

Fig. 4.

Verification of the negativity of the orbital derivative of the function V computed in Objective III (left). Areas where it is not negative are plotted in red. Not only is it not negative outside of the periodic orbit, but also in a strip within it. By increasing the spatial resolution the area where the orbital derivative fails to be negative inside of the orbit becomes smaller (right).

By using larger values for the parameters xres,yres,T,n, the area within the periodic orbit where the orbital derivative fails to be negative is effectively reduced to an arbitrary small neighbourhood of the equilibrium at the origin, cf. [12, Figs. 1,2]. Using a compiled programming language like C++ also makes these computations very fast.

4 Conclusion

We presented a new project for undergraduate students in Numerical Analysis. It is not difficult to solve the problem and thus the average student should be able to do so without too much difficulties. The solution, i.e. the function plotted, is somewhat surprising and intended to awaken the interest of the better students. The solution to the project is given and the results and the underlying theory are discussed in some detail. Matlab code is given for all objectives of the project together with code for a more sophisticated verification of the results. Further, numerous references to the underlying theory are given.

Let us compare the presented project, in terms of pedagogical value, to the larger projects from the textbook [20] of the course. The textbook contains, apart from numerous Exercises and Computer Problems, eleven so-called Reality Checks, which are larger project of comparable scale and complexity to the project presented. These Reality-Checks are very well designed and include topics such as the kinematics of the Steward Platform frequently used in flight-simulators, positioning using GPS, motion control in computer-aided modeling, and a simple audio codec. These Reality Checks have been used as assignments for the larger projects in the course discussed in the paper with a great success. The difference, however, between those Reality Checks and the project described, is that the Reality Checks inspire the better students to become interested in the technology behind the project, e.g. flight simulators, GPS, audio codecs, but not in mathematics. The teacher of the course has had lively discussion with students about these technological topics during and after the Reality Checks projects, but rarely about mathematics. A possible drawback of the proposed project is that the theory behind the project has little direct connection to technology. While a project on, say the Stewart Platform, might invoke the interest of an engineering student in numerical root-finding of nonlinear equations because she/he finds the steering of a flight simulator fascinating, she/he will just routinely solve the project proposed in this paper without gaining any interest in numerical analysis at all.

The advantage of the proposed project is that it inspires some of the mathematics students to get interested in the theory of dynamical systems, a topic that they had only known indirectly from studying differential equations. Further, as Lyapunov functions are a mathematical extension of the concept of dissipative energy from physics, some of the physics students were particularly interested as well. This allowed to the teacher to discuss attractors, repellers, chaos, and complete Lyapunov functions with interested students.

Appendix

Below is an implementation in Matlab to a solution of the project.

Program I:

Program II:

Program III: Modification of TraLengths to additionally verify the validity of the conditions for a Lyapunov function (decrease of the orbital derivative).

Add the following code to the function TraLengths between lines 19 and 20.

Further, add the following functions after the code for the function TraLengths