R TUTORIAL for the First Course. Part 3: Numerical Methods

Introduction

It would be very nice if discrete models provide calculated solutions to differential (ordinary and partial) equations exactly, but of course they do not. In fact in general they could not, even in principle, since the solution depends on an infinite amount of initial data. Instead, the best we can hope for is that the errors introduced by discretization will be small when those initial data are resonably well-behaved.

In this chapter, we discuss some simple numerical method applicable to first order ordinary differential equations in normal form subject to the prescribed initial condition:

\[ y' = f(x,y), \qquad y(x_0 ) = y_0 . \]

We always assume that the slope function \( f(x,y) \) is smooth enough so that the given initial value problem has a unique solutions and the slope function is differentiable so that all formulas make sense. Since it is an introductory, we do not discuss discretization errors in great detail---it is a subject of another course. A user should be aware that the magnitude of global errors may depend on the behavior of local errors in ways that ordinary analysis of discretization and rounding errors cannot predict. The most part of this chapter is devoted to some finite difference methods that proved their robastness and usefullness, and can guide a user how numerical methods work. Also, presented codes have limited applications because we mostly deal with uniform truncation meshes/grids that practically are not efficient. We do not discuss adaptive Runge--Kutta--Fehlberg algorithms of higher order that are used now in most applications, as well as spline methods. Instead, we present some series approximations, incuding Adomian Decomposition Method (ADM), as an introduction to application of recurrences in practical life.

The finite diference approximations have a more complecated "physics" than the equations they are designed to simulate. In general, finite recurrences usually have more propertities than their continuous anologous counterparts. However, this arony is no paradox for finite differences are used not because the numbers they generate have simple properties, but because those numbers are simple to compute.

Solving ODEs Numerically
The main R package deSolve to solve initial value problems (IVP) of:

The implementation includes stiff and nonstiff integration routines based on the ODE-PACK FORTRAN codes. Package deSolve contains several IVP ordinary differential equation solvers, that belong to the most important classes of solvers. It also has a variety of Runge--Kutta methods for solving stiff and not stiff equations.

Note that another package, bvpSolve provides methods to solve boundary value problems.

 

  1. Schiesser, W.E., An Introductory Comparison of Matlab and R, Lehigh University, Bethlehem, USA.
  2. Soetaert, K, Cash, J., Mazzia, F., Solving Differential Equations in R, 2012, Springer, Heidelberg New York Dordrecht London