Vector Fields and Solutions to

Ordinary Differential Equations using Octave

Andreas Stahel

16th December 2009

Contents
1 Vector fields 1
1.1 Vector field for the logistic equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Solutions of ordinary differential equations with lsode() . . . . . . . . . . . . . . . . . . 2
1.3 Vector field for the equation of a damped pendulum . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Solution to the pendulum equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 The ODE Package on Octave-Forge 6

2.1 Examples for ode23(), ode45() and ode78() . . . . . . . . . . . . . . . . . . . . . . 7

3 Codes from Lecture Notes 8

3.1 Using RK45(), Runge-Kutta adaptiv . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.2 Using ode Runge(), Runge-Kutta with fixed step . . . . . . . . . . . . . . . . . . . . . . 9
3.3 Using lsode() . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

In these notes you find the methods to generate vector fields with the help of Octave . In addition one of
the methods to solve ordinary differential equations is illustrated.
The notes assume that you have a recent version of Octave installed.

1 Vector fields
A vector field in the plane R2 is determined by a vector function
!
F1 (x1 , x2 )
F~ (~x) =
F2 (x1 , x2 )

At each point ~x ∈ R2 the vector F~ (~x) is attached. One of the possible applications of vectors fields is the
visualization of solution of ordinary differential equations.

1.1 Vector field for the logistic equation

The vector field for the logistic differential equation
d
x(t) = x(t) − x2 (t)
dt

1
1 VECTOR FIELDS 2

is given by !
1
F~ (t, x) =
x − x2
To visualize this field we have to

• choose the domain to be used

• define the functions for the vector field

• choose how many vectors to draw

• compute the vectors at the chosen points

• rescale the length of the vectors

• generate the graphic on screen

• save the result in a file

The above tasks can be performed with Octave , as shown in the code below with the resulting Figure 1.
Octave
t = 0 : 0 . 5 : 8 ; %% choose t h e time v a l u e s
x = 0 : 0 . 2 : 2 ; %% choose t h e x−v a l u e s

%% d e f i n e t h e f u n c t i o n
f u n c t i o n dx = l o g i s t i c ( x , t )
dx = x−x ˆ 2 ;
endfunction

%% c r e a t e t h e zero v e c t o r f i e l d
n t = l e n g t h ( t ) ; nx = l e n g t h ( x ) ;
V1 = z e r o s ( nt , nx ) ; V2 = z e r o s ( nt , nx ) ;

%% compute t h e v e c t o r f i e l d
for i = 1: nt
f o r j = 1 : nx
V1( i , j ) = 1 ; V2( i , j ) = l o g i s t i c ( x ( j ) , t ( i ) ) ;
endfor
endfor

%% choose t h e s c a l e f a c t o r
scalefactor = 0.2;

figure (1);
q u i v e r ( t , x , V1’ ,V2’ , s c a l e f a c t o r )
axis ([0 ,8 ,0 ,2]);
x l a b e l ( ’ time t ’ ) ; y l a b e l ( ’ p o p u l a t i o n x ’ )
g r i d on
%% p r i n t ( ” L o g i s t i c . png ” )
%% p r i n t ( ” L o g i s t i c . eps ” )

1.2 Solutions of ordinary differential equations with lsode()

Octave provides a selection of commands to solve ordinary differential equations, including systems. Here
we only use lsode() to generate solutions for the differential equations arising from the above vector
fields.

SHA 16-12-09
1 VECTOR FIELDS 3

1.5

population x 1

0.5

0
0 1 2 3 4 5 6 7 8
time t

Figure 1: Vector field for the logistic equation

To solve the initial value problem

d
x(t) = x(t) − x2 (t) with x(0) = 0.2
dt
we proceed as follows.

• Choose the values of time t for which the solution is to be computed.

• Compute the solution with the command lsode().

• We may compute more solution by using different initial conditions, e.g. x(0) = 0.4 and x(0) = 2.0.

• Plot the solutions.

• The solutions and the vector field can be shown in one graphic.

Find the result in Figure 2.

Octave
T = 0:0.1:8;

X1 = l s o d e ( ’ l o g i s t i c ’ , 0 . 2 , T ) ;
X2 = l s o d e ( ’ l o g i s t i c ’ , 0 . 4 , T ) ;
X3 = l s o d e ( ’ l o g i s t i c ’ , 2 . 0 , T ) ;

figure (3);
hold o f f
p l o t (T , [ X1, X2, X3 ] ) ;
g r i d on ; a x i s ( [ 0 , 8 , 0 , 2 ] ) ; hold on
q u i v e r ( t , x , V1’ ,V2’ , 0 . 4 )

1.3 Vector field for the equation of a damped pendulum

For atonomous differential equations the first component of the vector field equals 1. For general problems
this is not the case. As an example we convert the differential equation for a damped pemdulum

ÿ(t) = −k y(t) − α ẏ(t)

SHA 16-12-09
1 VECTOR FIELDS 4

1.5

0.5

0
0 1 2 3 4 5 6 7 8

Figure 2: Solutions and vector field for the logistic equation

into a system of diffeential equations of order one

! !
d y(t) v(t)
=
dt v(t) −k y(t) − α v(t)

and examine the resulting vector field. With the numerical values k = 1 and α = 0.1 we have to examine
the vector field !
1
F~ (y, v) =
−y − 0.1 v
The code to be used is very similar to the logistic equations and thus shown without any further explanation.
Octave
y = −1:0.1:1; %% choose t h e y v a l u e s
v = −1:0.1:1; %% choose t h e v−v a l u e s

%% d e f i n e t h e f u n c t i o n
f u n c t i o n dy = F1 ( y , v )
dy = v ;
endfunction

f u n c t i o n dv = F2 ( y , v )
k = 1 ; alpha = 0 . 1 ;
dv = −k∗y−alpha ∗v ;
endfunction

%% c r e a t e t h e zero v e c t o r f i e l d
ny = l e n g t h ( y ) ; nv = l e n g t h ( v ) ;
V1 = z e r o s ( ny , nv ) ; V2 = z e r o s ( ny , nv ) ;

%% compute t h e v e c t o r f i e l d
f o r i = 1 : ny
f o r j = 1 : nv
V1( i , j ) = F1 ( y ( i ) , v ( j ) ) ; V2( i , j ) = F2 ( y ( i ) , v ( j ) ) ;
endfor
endfor

SHA 16-12-09
1 VECTOR FIELDS 5

%% choose t h e s c a l e f a c t o r
scalefactor = 0.5;

figure (1);
q u i v e r ( y , v , V1’ , V2’ , s c a l e f a c t o r )
g r i d on ; a x i s ( [ min ( y ) , max( y ) , min ( v ) , max( v ) ] )
x l a b e l ( ’ p o s i t i o n x ’ ) ; y l a b e l ( ’ speed v ’ )

The only additional feature is the explicit choice of the domain shown in the graphic by the command
axis().
As can be seen in the left of Figure 3 the vectors at the origin have a small length and thus it is difficult to
read of the direction. This can be improved by normalizing all the vectors to have length 1 and then rescale
them to have a good length for visalisation. This is done in the code below witrh the result on the right in
Figure 3.
Octave
figure (2);
Vlength = s q r t (V1. ˆ 2 +V2 . ˆ 2 ) ;
V1n = V1 . / Vlength ; V2n = V2 . / Vlength ;
q u i v e r ( y , v , V1n ’ , V2n ’ , s c a l e f a c t o r )
g r i d on ; x l a b e l ( ’ p o s i t i o n x ’ ) ; y l a b e l ( ’ speed v ’ )
a x i s ( [ min ( y ) , max( y ) , min ( v ) , max( v ) ] )

1 1

0.5 0.5
speed v

speed v

0 0

-0.5 -0.5

-1 -1
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
position x position x

Figure 3: Standard and normalized vector field for the pendulum equation

1.4 Solution to the pendulum equation

To determine solutions of the system of differential equations
! ! ! !
d y(y) v(t) y(0) 0
= with =
dt v(t) −y(t) − 0.1 v(t) v(0) 0.9

we have to rewrite the function to be used by lsode().

Octave
T = 0:0.05:25;

f u n c t i o n dy = ODEPend( y , t )

SHA 16-12-09
2 THE ODE PACKAGE ON OCTAVE-FORGE 6

k = 1 ; alpha = 0 . 1 ;
dy = [ y ( 2 ) ; −k∗y(1)− alpha ∗y ( 2 ) ] ;
endfunction

YV = l s o d e ( ’ODEPend ’ , [ 0 ; 0 . 9 ] , T ) ;

figure (3);
p l o t (T ,YV) ;
g r i d on
x l a b e l ( ’ time ’ ) ; y l a b e l ( ’ p o s i t i o n and v e l o c i t y ’ ) ;

figure (4)
hold o f f
q u i v e r ( y , v , V1n ’ , V2n ’ , s c a l e f a c t o r )
g r i d on
a x i s ( [ min ( y ) , max( y ) , min ( v ) , max( v ) ] )
hold on
p l o t (YV( : , 1 ) ,YV( : , 2 ) , ’ r ’ )
xlabel ( ’ position ’ ) ; ylabel ( ’ velocity ’ ) ;
hold o f f

1
1

0.5 0.5
position and velocity

velocity

0 0

-0.5 -0.5

-1 -1
0 5 10 15 20 25 -1 -0.5 0 0.5 1
time position

Figure 4: One solution and the vector field for the pendulum equation

2 The ODE Package on Octave-Forge

The standard command in Octave to solve ordinary differential equations is lsode(), as used in the
previous section. On OctaveForge at [Link] you can find an package
with additional commands to solve ODEs.

• Dowload the package odepkg in a local directory

• Launch Octave and install the package with the command

pkg install [Link]

• Depending on your local configuration you might have to load the package with
pkg load odepkg

SHA 16-12-09
2 THE ODE PACKAGE ON OCTAVE-FORGE 7

• With this package many commands from MATLAB are now available in Octave , and some additional
commands. Table 1 shows some of the commands.

Solver Description
lsode() efficient, addaptive solver based on Hindmarsh’s work ([Hind93])
ode23() addaptive, explicit,solver, based on Heun’s method
ode45() addaptive, explicit solver, based on Runge-Kutta method of order 4, resp. 5
ode54() Runge Kutta solver
ode78() addaptive, explicit solver, based on Runge-Kutta method of order 7, resp. 8
ode2r() solver for stiff problems, based on Hairer and Wanner’s code
ode5r() solver for stiff problems, based on Hairer and Wanner’s code

Table 1: Octave commands to solve ordinary differential equations

The package provides many additional commands, some of them shown in Table 2. Find a description
of the algorithms, their advantages and disadvantages in the file [Link], to be found in the sub-
directory with the package.

Command Description
odeexamples() launch demos for ordinary differential equations
ode23d() solver for delay differential equations
ode45d() solver for delay differential equations
ode78d() solver for delay differential equations

Table 2: Additional commands in the ODE package

2.1 Examples for ode23(), ode45() and ode78()

As an example we consider again the pendulum equation
! ! ! !
d y(y) v(t) y(0) 0
= with =
dt v(t) −y(t) − 0.1 v(t) v(0) 0.9
Thus a function file with this function may be generated and then used by all examples below.
ODEPend.m
f u n c t i o n dy = ODEPend( t , y )
k = 1 ; alpha = 0 . 1 ;
dy = [ y (2); −k∗y(1)− alpha ∗y ( 2 ) ] ;
endfunction

As a general rule we first choose the parameters for the solver. As a first example we choose a relative
tolerance of 10−3 and an absolute tolerance of 10−3 . We want a graph to be generated while the differential
equation is solved.
Octave
vopt = o d e s e t ( ” RelTol ” , 1e−3, ”AbsTol ” , 1e −3,”NormControl ” , ”on ” ,\
” OutputFcn ” , @odeplot ) ;

Then we use ode78() to solve the system of equations.

Octave

SHA 16-12-09
3 CODES FROM LECTURE NOTES 8

ode78 (@ODEPend, [0 2 5 ] , [0 0 . 9 ] , vopt ) ;

The resulting animation and final solution will look rather ragged. Since the Runge–Kutta algorithm of
order 7 is very efficient, only very few points have to be computed. Thus the graphic does not look nice, but
the numerical results are reliable. The command
Octave
ode23 (@ODEPend, [0 2 5 ] , [0 0 . 9 ] , vopt ) ;

will generate a nice looking graph, but require more computation time.

With all of the above codes the numerical values will not be returned and are thus not available for
further computation. Use the code below if further computations have to be performed. You will find a plot
of position and velocity as function of time and a phase plot.
Octave
vopt = o d e s e t ( ” RelTol ” , 1e−10, ”AbsTol ” , 1e −10,”NormControl ” , ”on ” ) ;
[ t , y ] = ode78 (@ODEPend, [0 2 5 ] , [0 0 . 9 ] , vopt ) ;
figure (1); plot ( t , y );
x l a b e l ( ’ time ’ ) ; y l a b e l ( ’ p o s i t i o n and v e l o c i t y ’ ) ; g r i d on
figure (2); plot (y (: ,1) , y ( : , 2 ) ) ;
x l a b e l ( ’ p o s i t i o n ’ ) ; y l a b e l ( ’ v e l o c i t y ’ ) ; g r i d on

3 Codes from Lecture Notes

To illustrate the codes presented in the lecture notes we use the example of a diode circuit examined in the
notes. The behavior of the diode is given by a function
(
0 for u ≥ −us
i = D (u) =
RD (u + us ) for u < −us

Based on Kirchhoff’s law we find the differential equations

1
u̇h = (−D (uh − uin ) + D (uout − uh ))
C1
1
u̇out = u̇in − D (uout − uh )
C2
We define the initial conditions and the two function in a script file. We examine an input voltage of
uin (t) = 10 cos(t).
Octave
Tend = 30; u0 = [ 0 ; 0 ] ;

f u n c t i o n c u r r = Diode ( u )
Rd = 10; us = 0 . 7 ;
i f ( u>=−us ) c u r r =0;
e l s e c u r r =Rd∗( u+us ) ;
endif
endfunction

function y = c i r c u i t (u , t )
C1 = 1 ; C2 = 1 ;
y = [ −1/C1∗( Diode ( u(1)−10∗ s i n ( t ))−Diode ( u(2)−u ( 1 ) ) ) ;
10∗ cos ( t ) −1/C2∗Diode ( u(2)−u ( 1 ) ) ] ;
endfunction

SHA 16-12-09
3 CODES FROM LECTURE NOTES 9

3.1 Using RK45(), Runge-Kutta adaptiv

We can use the adaptive Runge-Kutta algorithm with relative and absolute tolerance of 10−5 and generate
the plot by
Octave
t 0 = cputime ( ) ;
[ t , u ] = rk45 ( ’ c i r c u i t ’ , 0 , Tend , u0 , 1 e−5,1e −5); % Runge Kutta a d a p t i v
t i m e r = cputime ()− t 0
figure (1);
plot ( t , u (: ,2) ’. ’)
g r i d on ; x l a b e l ( ’ time ’ ) ; y l a b e l ( ’ t e n s i o n ’ ) ;

The result in Figure 5 seems reasonable, but it took 8 seconds of CPU time to compute.

3.2 Using ode Runge(), Runge-Kutta with fixed step

We can compare the result with a Runge-Kutta calculation with 100 steps, resulting in Figure 5.
Octave
t F i x = l i n s p a c e ( 0 , Tend , 1 0 0 ) ;
t 0 = cputime ( ) ;
[ tFix , uFix ] = ode Runge ( ’ c i r c u i t ’ , tFix , u0 , 1 ) ; % Runge Kutta
t i m e r = cputime ()− t 0
p l o t ( tFix , uFix ( : , 2 ) , t , u ( : , 2 ) )
g r i d on ; x l a b e l ( ’ time ’ ) ; y l a b e l ( ’ t e n s i o n ’ ) ;
legend ( ’ u f i x ’ , ’ u adapt ’ )

It took only 0.15 seconds, but the solution is ragged at the turning points. This can be improved by using 10
intermediate steps between the output times. Use
Octave
[ tFix , uFix ] = ode Runge ( ’ c i r c u i t ’ , tFix , u0 , 1 0 ) ; % Runge Kutta

to find a competitive solution in 1.3 seconds. This is one of the rare occasion where a fixed size algorithm
outperforms an adaptive algorithm.

3.3 Using lsode()

We can compare the above result with the performance of lsode(). Unfortunately the arguments t and
u have to be given in reverse order for lsode() and the above codes and thus the header of the function
circuit has to be modified slightly. We have to specify the tolerances. With a computation time of 0.48
seconds we find a good solution. This illustrates the quality of the algorithm in lsode.m.
DoubleTensionLSODE.m
Tend = 30; u0 = [ 0 ; 0 ] ;

f u n c t i o n c u r r = Diode ( u )
Rd = 10; us = 0 . 7 ;
i f ( u>=−us ) c u r r =0; e l s e c u r r =Rd∗( u+us ) ; endif
endfunction

function y = c i r c u i t (u , t )
C1 = 1 ; C2 = 1 ;
y = [ −1/C1∗( Diode ( u(1)−10∗ s i n ( t ))−Diode ( u(2)−u ( 1 ) ) ) ;
10∗ cos ( t ) −1/C2∗Diode ( u(2)−u ( 1 ) ) ] ;
endfunction

t = l i n s p a c e ( 0 , Tend , 1 0 0 ) ;

SHA 16-12-09
REFERENCES 10

30
u fix
u adapt
25

20

15
tension

10

-5
0 5 10 15 20 25 30
time

Figure 5: Solution for the diode circuit with Runge-Kutta, fixed step and adaptive

l s o d e o p t i o n s ( ” a b s o l u t e t o l e r a n c e ” ,1 e −5);
l s o d e o p t i o n s ( ” r e l a t i v e t o l e r a n c e ” ,1 e −5);

t 0 = cputime ( ) ;
u = l s o d e ( ’ c i r c u i t ’ , u0 , t ) ;
t i m e r = cputime ()− t 0
plot ( t , u (: ,2))
g r i d on ; x l a b e l ( ’ time ’ ) ; y l a b e l ( ’ t e n s i o n ’ ) ;

In the previous section the codes and data files in Table 3 were used.

filename function
ode Euler.m algorithm of Euler, fixed step size
ode Heun.m algorithm of Heun, fixed step size
ode Runge.m algorithm of Runge-Kutta, fixed step size
rk45.m adaptive Runge Kutta algorithm
DoubleTension.m sample code for the diode circuit
DoubleTensionLSODE.m using lsode()

Table 3: Codes and data files

References
[Hind93] A. C. Hindmarsh and K. Radhakrishnan. Description and Use of LSODE, the Livermore Solver
for Ordinary Differential Equations. NASA, 1993.

SHA 16-12-09