Tutorial 9: numeric solutions of ODEs
Objectives:
- Use functions and the deSolve package to calculate numeric solutions of ODEs
- Use functions to graph the defining functions of ODEs
- Pass function names to other functions (optional)
Numeric solution of differential equations
We will use the package deSolve to calculate numeric solutions of ODEs. First, we need to create the R function that defines the derivative in the differential equation, in other words, the function \(f(x,t)\) in the generic first-order ODE
\[ \frac{dx}{dt} = f(x,t)\]
This R function must have three inputs: t, x, parms, in that order, representing the time variable, the dependent variable x, and the vector of parameters parms. Let us illustrate this on the linear population model with birth rate b and death rate d:
\[ \frac{dN}{dt} = bN-dN \]
To solve this ODE with define the function that calculates the function on the right-hand-side (bN-dN). Notice that even though the function does not depend on time, t still must be the first input argument of the function:
The first two lines extract the two parameters b and d from the vector of parameter values, the next line calculate the value of the derivative, and the last one returns it.
Third, we assign the parameter values to the vector, create a time vector on which to solve the ODE, and assign the initial condition(s):
Finally, we are ready to call the function ode that will do all the work to solve this ODE using the function pop_funk:
The output is a data frame, which means that you can use the data= option in plot to make your graph of the solution as a function of time:
This solution, like all numeric solutions of ODEs, is an approximation of the exact (analytic) solution. In this case, we can find (and verify) the exact solution of this model:
\[ N(t) = N(0)e^{(b-d)t}\]
Thus we can plot the analytic solution along with the numeric solution to see how far off they are:
You can see that the ODE solver is so clever that the numeric solution appears identical to the exact solution, even though there is always some degree of error in numeric solutions of ODEs.
Plotting defining functions of ODEs
Suppose we modify the defining function of the ODE to include a constant term C:
\[ \frac{dN}{dt} = bN-dN + C\]
The R function for the ODE can be defined as follows:
To analyze the ODE graphically, let us create a plot of the defining function \(f(N)\) over a range of values of \(N\). This requires choosing a range that includes all the zeros of the function \(f(N)\), which are the fixed points of the ODE. So if we let \(b=0.3\), \(d=0.32\), and \(C=1\), the function \(f(N) = -0.02N+1\) has a zero at 50, so we can assign the range of values of N from 0 to 100 and make a graph of the function over this range:
Note that we need to use the function unlist() to turn dNdt into a regular vector from a list (the list structure is necessary for it to work with the function ode). Also note that we had to define the vector time even though it is not used in the calculation because it is an input of the function pop_funk2m again because itβs required by ode.
The plot of the defining function shows the rate of change of the solution (dN/dt) as a function of N. For population values below the fixed point of 50, solutions grow, while for N>50 solutions decay, both converging to the asymptotic value of 50. This can be shown by plotting several solutions obtained by calling ode:
Calling functions using strings (optional)
For the curious, here is a way to specify and call a function based on a given character string. You can see that calling the new_fun()is the same as calling the original function blah():
This is very useful if you want to pass the name of a function as a string (e.g. blah) to another function (e.g. my_funk), so then it can be used to call the specified function from within my_funk. This allows you to write a general function that can call any number of functions and perform the same calculations with them.