Numerical Methods for Differential Equations
Forest W. Arnold
November 2018
1
Typeset in LATEX.
Copyright © 2018 Forest W. Arnold
This work is licensed under the Creative Commons Attribution-Noncommercial-Share
Alike 3.0 United States License. To view a copy of this license, visit
[Link] or send a letter to
Creative Commons, 171 Second Street, Suite 300, San Francisco, California, 94105,
USA. You can use, print, duplicate, share this work as much as you want. You can base
your own work on it and reuse parts if you keep the license the same.
Trademarks
TI-Nspire is a registered trademark of Texas Instruments, Inc.
2
1 Purpose
The objectives of this article are to
1. describe and demonstrate TI-Nspire’s built-in numerical methods for solving
initial-value problems;
2. describe and implement the Runge-Kutta 4th order numerical method (RK4);
3. describe and implement the Runge-Kutta-Fehlberg 4th order adaptive numerical
method (RK45);
4. compare the built-in RK23 method, the RK4 method, and the RK45 method; and
5. show how to graph numerical solutions with TI-Nspire CAS’s interpolate()
function.
2 Description
Although there are a wide variety of analytical techniques for solving differential equa-
tions, most real-world differential equations can not be solved with those techniques.
Instead, numerical methods for finding approximate solutions must be used. Both the
CAS and non-CAS versions of TI-Nspire have two built-in numerical functions for
solving differential equations: euler() and rk23(). As implied by the names, the
euler() function implements Euler’s method and the rk23() function implements
the Runge-Kutta second-order method with third-order error control. Both of these
methods accept either a single first-order differential equation or a system of differen-
tial equations as input and return a matrix of values as output.
Both Euler’s method, its derivative methods, and the Runge-Kutta methods work by
creating a mesh of grid points and approximating values using estimates of the slope
between successive grid points, then using the estimated slope to approximate the co-
ordinates of the succeeding grid point. The main difference in these methods is how
the slope between grid points is estimated. For the Euler method, the slope between
the grid points (xi , yi ) and (xi+1 , yi+1 ) is estimated by simply evaluating f (xi , yi ), where
f (x, y) is the formula for the differential equation. For the Runge-Kutta methods, the
slope is approximated using a weighted average of estimates of the slopes at several
points in the interval [xi , xi+1 ]. The order of a particular Runge-Kutta method indicates
the number of points used to estimate the slope between the endpoints of the interval.
For example, a fourth-order Runge-Kutta method estimates the slope using a weighted
average of four slope estimates.
Euler’s method is the simplest method for approximating solutions to initial-value dif-
ferential equations. Because of its simplicity, it can be used to illustrate the basic tech-
nique and the steps involved in finding approximate solutions. Figure 1 graphically
illustrates how Euler’s method approximates values at each grid point.
3
Figure 1: Euler’s Method
Given an initial value x, y, a solution domain [a, b], and a step size h (the length of
the interval between values of the independent variable), Euler’s method starts with the
specified value of the equation y(x0 ) = y0 . It then calculates succeeding values in steps:
yi+1 = yi + slope · h where slope is the derivative of the equation at the point (xi , yi )
and h is the step size. The first three steps of this process are shown in Figure 1. The
blue curve is the graph of the actual (unknown) solution, the red lines are the tangent
lines at the grid points (xi , yi ), and the red dots are the linear approximations at each
grid point (xi+1 , yi+1 ).
Since linear approximations are used to find values of successive points, unless the
actual solution is a linear function, except for the initial value the calculated values will
differ from the actual unknown values. The differences are the approximation errors,
and the overall error accumulates with each step. For methods with fixed step sizes
such as Euler’s method or the RK4 method, the only way to reduce the magnitude of
the errors is by decreasing the step size. Additionally, there is no straight-forward way
to determine the step size required to achieve a specified error tolerance.
Adaptive methods such as rk23() and rk45() (implemented below) attempt to satisfy
a specified tolerance value by automatically modifying the step size. If the specified
tolerance at a step is not met, the step size is reduced and the calculations for the step
are repeated using the modified step size. These methods can also increase step sizes,
thus reducing the total number of steps needed.
3 Comparing euler() and rk23()
This section demonstrates using TI-Nspire’s built-in functions for numerically approx-
imating the solution to a first-order initial-value problem and compares the results of
the two methods. The equation to approximate is y0 = y + e2x with initial value y(0) = 1
over the interval [0, 1] using a step size of 0.1. The solution of this equation is y = e2x .
4
3.1 Approximating the Solution with euler()
The TI-Nspire function euler() by default simply uses the specified step size, varstep,
to determine how many steps to take to find the overall approximation. This function
also accepts an optional argument, eulerstep, which specifies how many steps to take
between each varstep.
The following command executes euler() with varstep equal to 0.1 without the op-
tional eulerstep argument:
The returned value is a matrix, with row 1 entries equal to the x values and row 2
entries equal to the y values. To graph the points, the x and y entries are extracted as
lists. To graph the curve through the points using the interpolate() function, a list
of the slopes (derivatives) at each point is created using the differential equation and
the constraint operator, | to substitute x and y values in the equation:
To show the effect of increasing the number of internal steps, the following command
executes euler() with varstep equal to 0.1 with the optional eulerstep argument equal
to 2 (twice as many steps):
The lists of x and y values are extracted from the matrix and a list of the slopes at
each point is created:
5
Figure 2 graphically displays the results of the two approximations and compares the
graphs of the approximate solutions with the graph of the actual solution.
Figure 2: Approximations Using the euler() Function
It is clear from the graph in Figure 2 that Euler’s method does not provide very precise
approximations unless a large number of steps are used. TI-Nspire’s Lists & Spread-
sheet application can be used to calculate the absolute error of approximation, which
equals the absolute value of the difference between the actual value and the approxi-
mate value. Figure 3 shows a spreadsheet with the calculated errors.
6
Figure 3: Approximation Errors with the euler() Function
3.2 Approximating the Solution with rk23()
The TI-Nspire function rk23() by default uses a tolerance value of 0.001 for the error
tolerance at each step. The default tolerance value can be over-ridden with an optional
argument, diftol, to specify a different tolerance value.
The following command executes rk23() with a step size equal to 0.1 with the de-
fault tolerance value, 0.001:
Next, the rows of x and y values are extracted from the matrix as lists for plotting
using the Graph application’s scatter plot functionality:
To show the result of using the optional diftol value, the function is executed with a
diftol value of 1.e−5 :
7
Figure 4 displays the graph of the actual solution and the approximate solutions from
executing rk23() with both the default tolerance value and a tolerance value of 1.e−5 :
Figure 4: Approximations Using the rk23() Function
As Figure 4 shows, the results from executing rk23() with both the default tolerance
value and a stricter tolerance value approximate the actual solution fairly well. Figure
5 displays the absolute error of approximation for the rk23() function.
8
Figure 5: Approximation Errors with the rk23() Function
These approximations are also much better than the approximations from the euler()
function. Figure 6 shows a spreadsheet comparing the absolute errors from the rk23()
and euler() functions.
Figure 6: Comparison of Approximation Errors
4 Higher Order Runge-Kutta Methods
Even though the Runge-Kutta methods perform multiple function calls per step, these
methods provide better performance and accuracy with larger step sizes than methods
such as Euler’s method or Modified Euler’s method. For this reason, the most-used nu-
merical methods for approximating non-stiff first-order ordinary differential equations
are the order 4 Runge-Kutta method and the adaptive order 4/order 5 Runge-Kutta-
Fehlberg method. Neither of these methods are provided with TI-Nspire, but both can
be easily implemented with TI-Nspire’s Program Editor.
9
4.1 The Order 4 Runge-Kutta Method
4.1.1 Description of the Method
The order 4 Runge-Kutta method approximates successive points (xi+1 , yi+1 ) at each
step using the average of four slope estimates per step interval. Given a step size h and
xi+1 = xi + ih where xi and xi+1 are the endpoints of the interval for step i, the formula
for calculating wi+1 , the approximation for yi+1 is
k1 = h f (xi , wi )
h k1
k2 = h f (xi + , wi + )
2 2
h k2
k3 = h f (xi + , wi + )
2 2
k4 = h f (xi + h, wi + k3)
1
wi+1 = wi + (k1 + 2 · k2 + 2 · k3 + k4)
6
4.1.2 Algorithm
The alogorithm 1 for the order 4 Runge-Kutta method is:
Runge-Kutta Order 4 Algorithm
Given a reference to a continuous function f (x, y), endpoints of the inter-
val of integration a and b, an initial value y(a) = al pha, and n, the number
of steps, determine an approximate solution of the initial value problem
y0 = f (x, y), a ≤ x ≤ b, y(a) = al pha.
create 2 by n+1 matrix of return values
set matrix[1, 1] = a
set matrix[2, 1] = al pha
set h = b−a
n
for i = 2, i <= n + 1
set im1 = i − 1
set xm1 = matrix[1, im1]
set wm1 = matrix[2, im1]
set k1 = h · f (xm1, wm1)
set k2 = h · f (xm1 + 2h , wm1 + k1
2)
set k3 = h · f (xm1 + h2 , wm1 + k2
2)
set k4 = h · f (xm1 + h, wm1 + k3)
set matrix[1, i] = xm1 + h
set matrix[2, i] = wm1 + 16 (k1 + 2 · k2 + 2 · k3 + k4)
end
return matrix
1 Algorithm 5.2, Burden and Faires, pp 288-289
10
4.1.3 Programming the Method
Writing the code for the method with TI-Nspire’s programming language is almost a
direct translation of the pseudo-code algorithm for the method. To code the method:
Step 1: Add a Program Editor page to the document by selecting the ”Insert - Pro-
gram Editor - New” menu item,
Step 2: In the pop-up dialog, enter rk4 for the Name entry, select Function for the
Type entry, and None for the Library Access entry.
Step 3: Type the program statements shown below in the program page. Be sure to
select the Check Syntax and Store after entering the program statements.
Step 4: Test the newly-implemented function by running it in a calculator page.
The listing for the code follows.
Define rk4(fname,a,b,alpha,numsteps)=
Func
© rk4(fname,a,b,alpha,numsteps) - numerically solve de using rk 4 algorithm
© input arguments:
© fname - the (string) name of the de of the form y’=f(t,y)
© a - left endpoint of the interval [a,b] for the solution
© b - right endpoint of the interval [a,b] for the solution
© alpha - initial value, y(a)=alpha
© numsteps - number of equally-spaced steps
© returns:
© a 2 by numsteps matrix with values of x in row 1 and y in row 2
©----------------------------------------------------
© This function is based on algorithm 5.2 on pages 288-289 of Burden
© and Faires "Numerical Analysis"
© Author: Forest W. Arnold
©====================================================
Local mat,h,i,im1,xm1,wm1,err,mode
Local hd2,k1,k2,k3,k4
If a>=b Then
Disp "Error: argument ’a’ must be less than argument ’b’!"
Return {}
EndIf
If numsteps<2 Then
Disp "Error: argument ’numsteps’ must be at least 2!"
Return {}
EndIf
mode:=setMode(5,2)
mat:=newMat(2,numsteps+1)
mat[1,1]:=a © xstart
mat[2,1]:=alpha © initial value at x start
h:=((b-a)/(numsteps)) © step size
For i,2,numsteps+1
im1:=i-1
xm1:=mat[1,im1] © previous x
wm1:=mat[2,im1] © previous w value
hd2:=((h)/(2))
k1:=h*#fname(xm1,wm1)
k2:=h*#fname(xm1+hd2,wm1+((k1)/(2)))
11
k3:=h*#fname(xm1+hd2,wm1+((k2)/(2)))
k4:=h*#fname(xm1+h,wm1+k3)
mat[1,i]:=xm1+h
mat[2,i]:=wm1+((k1+2*k2+2*k3+k4)/(6))
EndFor
setMode(5,mode)
Return mat
EndFunc
4.1.4 Testing the rk4() Function
Test the new function by approximating the solution to the initial-value problem y0 −
y = 0, y(0) = 1, 0 <= 1 using 10 steps (for a step size of 0.1). The actual solution of
the differential equation y0 = y is y = ex .
Check the solution by comparing the returned approximate y values with values of
y = ex .
4.2 The Adaptive Order4/Order5 Runge-Kutta-Fehlberg Method
4.2.1 Description of the Method
The adaptive Runge-Kutta-Fehlberg Method, developed by the German mathemati-
cian Erwin Fehlberg, automatically adjusts the step size for each iteration to achieve
a specified tolerance value. This is done by calculating both order 4 and order 5 ap-
proximations in each step, then comparing the order 5 and order 4 approximations to
estimate the ”expected error” of approximation. Only one additional function evalu-
ation per step is needed to calculate the order 5 approximation. The estimate of the
expected error is the difference between the order 4 and order 5 approximations. If the
expected error is larger than the tolerance value, the step size is reduced and the step
is repeated with the reduced step size. If the expected error is less than the tolerance
value, the step size is slightly increased for the next step.
4.2.2 Algorithm
The alogorithm 2 for the adaptive Runge-Kutta order4/order 5 method is:
2 Algorithm 5.3, Burden and Faires, pp 297-298
12
Runge-Kutta Order 4/Order 5 Algorithm
Given a reference to a continuous function f (x, y), endpoints of the interval of
integration a and b, an initial value y(a) = al pha, a tolerance value tol, a mini-
mum step size hmin, and a maximum step size hmax, determine an approximate
solution of the initial value problem y0 = f (x, y), a ≤ x ≤ b, y(a) = al pha.
set t = a; set w = al pha; set h = hmax; set i = 1; set f lag = 1
create a list tvals to return t values and a list yvals to return w values
set tvals[1] = t; set yvals[1] = w
while f lag = 1
set k1 = h · f (t, w)
set k2 = h · f (t + 14 h, w + 41 k1)
set k3 = h · f (t + 38 h, w + 32 3
k1 + 329
k2)
12
set k4 = h · f (t + 13 h, w + 2197 k1 − 7200
1932 7296
2197 k2 + 2197 k3)
set k5 = h · f (t + h, w + 439 3680 845
216 k1 − 8k2 + 513 k3 − 4104 k4)
set k6 = h · f (t + 2 h, w − 27 k1 + 2k2 − 2565 k3 + 1859
1 8 3544 11
4104 k4 − 40 k5)
1 1 128 2197 1 2
set r = h · | 360 k1 − 4275 k3 − 75240 k4 + 50 k5 + 55 k6|
Note: r=1/h times difference between order 4, order 5 approximations
if r <= tol
Note: approximation meets tolerance value
set t = t + h
25
set w = w + 216 k1 + 1408 2197 1
2565 k3 + 4104 k4 − 5 k5
set tvals[i] = t
set yvals[i] = w
set i = i + 1
end if
Note: Calculate step delta and see if step needs to be changed
41
set delta = 0.84 · tol r
if delta ≤ 0.1 set h = 0.1 · h
else if delta ≥ 4 set h = 4 · h
else h = delta · h
if h > hmax set h = hmax
Check for end of iteration
if t ≥ b set f lag = 0
else if t + h > b set h = b − t
else if h < hmin display a warning and return
end while
return a matrix of x, y values
4.2.3 Programming the Method
Add a Program Editor page to the document, name the function rk45 with a type of
Function and library access of None, then enter the following code in the editor.
13
Define rk45(fname,a,b,alpha,tol,hmin,hmax)=
Func
© rk45(fname,a,b,alpha,tol,hmin,hmax) - approximate a diffeq
© input arguments:
© fname - the (string) name of the diff eq of the form y’=f(x,y)
© a - the left endpoint of the interval, [a,b]
© b - the right endpoint of the interval, [a,b]
© alpha - initial value: y(a)=alpha
© tol - tolerance value specifying precision of each value
© hmin - minimum step size
© hmax - maximum step size
© returns:
© a 2 by n matrix with values of x in row 1 and y in row 2
©---------------------------------------------
© This function is based on algorithm 5.3 on pages 296-300 of
© Burden and Faires "Numerical Analysis"
© Author: Forest W. Arnold
©=============================================
Local t,w,h,del,flag,steps,mode
Local k1,k2,k3,k4,k5,k6
Local r,i,tvals,yvals,ypvals
If a>=b Then
Disp "Error: argument ’a’ must be less than argument ’b’!"
Return {}
EndIf
If hmin>=hmax Then
Disp "Error: argument ’hmin’ must be less than argument ’hmax’!"
Return {}
EndIf
mode:=setMode(5,2)
t:=a
w:=alpha
h:=hmax
i:=2
flag:=1
tvals:={t}
yvals:={w}
steps:=0
While flag=1
steps:=steps+1
k1:=h*#fname(t,w)
k2:=h*#fname(t+(1/4)*h,w+(1/4)*k1)
k3:=h*#fname(t+(3/8)*h,w+(3/32)*k1+(9/32)*k2)
k4:=h*#fname(t+(12/13)*h,w+(1932/2197)*k1-(7200/2197)*k2+
(7296/2197)*k3)
k5:=h*#fname(t+h,w+(439/216)*k1-8*k2+(3680/513)*k3-
(845/4104)*k4)
k6:=h*#fname(t+(1/2)*h,w-(8/27)*k1+2*k2-(3544/2565)*k3+
(1859/4104)*k4-(11/40)*k5)
r:=(1/h)*abs((1/360)*k1-(128/4275)*k3-(2197/75240)*k4+
(1/50)*k5+(2/55)*k6)
If r<=tol Then
t:=t+h
w:=w+(25/216)*k1+(1408/2565)*k3+(2197/4104)*k4-(1/5)*k5
tvals[i]:=t
yvals[i]:=w
i:=i+1
EndIf
del:=0.84*((tol/r))^(1/4)
If del<=0.1 Then
h:=0.1*h
ElseIf del>=4 Then
h:=4*h
Else
h:=del*h
EndIf
If h>hmax Then
h:=hmax
EndIf
14
If t>=b Then
flag:=0
ElseIf t+h>b Then
h:=b-t
ElseIf h<hmin Then
flag:=0
Disp "Warning: step size less than minimum at step ",steps,"!"
setMode(5,mode)
Return {tvals,yvals}
EndIf
EndWhile
setMode(5,mode)
Return {tvals,yvals}
EndFunc
4.2.4 Testing the rk45() Function
Test the new function by approximating the solution to the initial-value problem y0 −
y = 0, y(0) = 1, 0 <= 1 with a tolerance value of 1.e− 5. The actual solution of the
differential equation y0 = y is y = ex .
Check the solution by comparing the returned approximate y values with values of
y = ex .
4.3 Comparing the Runge-Kutta Functions
The initial-value problem y0 = y + e2x , y(0) = 1, 0 <= x <= 1 is approximated using
rk23(), rk4(), and rk45() to compare the performance and precision of each of
the three Runge-Kutta functions.
First, the differential equation is defined:
Then the three functions are executed and the lists of x, y values are extracted from
the returned matrices:
15
A graph with the lists of x, y values plotted as scatter plots against the plot of the actual
solution is shown in Figrue 7.
Figure 7: Approximations from rk23(),rk4(), and rk45()
16
As Figure 7 shows, all three of the Runge-Kutta functions generate good approxi-
mations to the solution y = e2x . Setting the Display Digits for the graph window to
Float10 and using the Graph application’s Trace All functionality, the approximations
for each x value can be displayed in the graph.
The absolute approximation errors can calculated in a Lists and Spreadsheet appli-
cation. Figure 8 is a spreadsheet comparing the absolute errors for the rk23() and
rk4() approximations.
Figure 8: Absolute Errors from rk23(),rk4() Approximations
Figure 9 is a spreadsheet with the absolute errors for the rk45() approximation.
Figure 9: Absolute Errors from the rk45() Approximation
From the above spreadsheets, rk4() approximates the solution better than rk23() and
rk45() generates the best approximation with fewer steps and with a global error less
than the specified tolerance value of 1.e−5 .
5 Graphing with the interpolate() Function
Given a list of x values and a list of corresponding y values, TI-Nspire’s built-in
interpolate() function can fit a curve through the x, y pairs using cubic interpo-
17
lation. This function can be used to plot the curve and can also be used to approximate
x, y values that are not in the lists.
The arguments to the interpolate() function are:
x - the x value for which the y value is to be returned.
xlist - a list of x values.
ylist - a list of y values corresponding to each x value in xlist.
yplist - a list of slopes at each x, y pair from xlist and ylist.
As shown in previous sections, the Runge-Kutta functions return matrices with x values
in row 1 and the corresponding y values in row 2. The rows can be extracted from the
matrices as lists, which can be used as input arguments to the interpolate() func-
tion. However, the function also requires a list of the slopes. For a first-order ordinary
differential equation, the required slopes can be calculated directly from the differential
equation using TI-Nspire’s constraint operator, |, as follows for the lists returned from
rk23(), rk4(), and rk45 for the equation f (x, y) = y + e2x :
Figure 10: Plots with the interpolate() Function
18
Figure 10 displays the plots of the cubic interpolants for the three Runge-Kutta func-
tions. To view each plot, place the cursor over the curves and use the Tab key to cycle
through the curves.
The interpolate() function can also be ”wrapped” in a user-defined convenience
function that returns an approximate y value corresponding to an input x value:
6 Summary
The first part of this article discussed and demonstrated TI-Nspire’s built-in functions
for numerically approximating solutions of ordinary differential equations, euler()
and rk23(). The approximations with the functions were compared and the absolute
errors were calculated.
The second part discussed the most commonly used numerical methods for differential
equations, rk4(), the Runge-Kutta order 4 method and rk45(), the Runge-Kutta-
Fehlberg order 4/order 5 adaptive method. Implementations of both these methods
using TI-Nspire’s Program Editor were presented and demonstrated, along with cal-
culations of the absolute errors of approximation with the methods.
The last part of this article discussed and demonstrated using TI-Nspire’s interpolate()
function for graphing lists of x, y values and finding approximate values of y using in-
terpolation.
References
[1] Burden, Richard L. and Faires, J. Douglas, Numerical Analysis, 9th Ed., Brook-
s/Cole Cengage Learning, 2011, pp. 67-71
[2] Gilat, Amos and Subramaniam, Vish, Numerical Methods for Scientists and En-
gineers, John Wiley & Sons,Inc., 2014, pp. 405-417
19