Ordinary Differential Equations part 1

To solve ordinary differential equations numerically, Mathcalc8 have two functions called 'rkode' and 'euode'. The first function using the method of 'Runge-Kutta'(rk) and the second function use the method of 'Euler' (eu). eu's method is simple and direct, so it is faster than rk but with less accuracy. Since both functions use exactly the same input parameters, we would only demonstrate 'rkode' in this post.
A general form of ordinary differential equations that we are going to solve in vector format are given as follows :

dx/dt = F(t , x) subject to 

Initial condition xt=t0 = x0

In the above equations, the variables and the function F are written in bold letter which means that they are vectors.
Consider the following example

Solve the equation for t from 0 to 5 :

dy/dx = -2x + y

Initial condition : xt=0 = 1

To solve this equation, we first need to rename the variables. The variable for the denominator of the differential operator is renamed as 't0'. The variable for the nominator of the differential operator is renamed as 'x0'. The renamed problem is shown as follows :

Solve the equation for t from 0 to 5 :

dx0/dt0 = -2t0 + x0

Initial condition : x0t0=0 = 1

Input 1 : r=rkode('-2*t0+x0','1',0,5,0.5)
Ans. set r = 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 1, 1.351, 1.283, 0.521, -1.384, -5.172, -12.065, -24.076, -44.523, -78.878, -136.158

The result is a solution matrix consists of two lists. The first element is a list r(0) is the time sequence from 0 to 1 with interval 0.5 which is set by the last parameter of the function. The second element is a matrix r(1) consists of one list which represents the solution x0 which is from 1 up to -136.158.
To get a plot of the solution :

Input 2 : d=r(1)
 

Input 3 : plot('line','r(0):d(0)',0,5,0.5)
 

The plot is shown below :

Actually, there is a formula for the solution which is :
y = 2 - e^x + 2x We can check this with our numerical solution r(1) by using r(0) as x and calculate y :

Input 4 : 2-exp(r(0))+2*r(0)
Ans. 1, 1.351, 1.282, 0.518, -1.389, -5.182, -12.086, -24.115, -44.598, -79.017, -136.413

It seemed that an interval setting of 0.5 is enough in this case for a match.
Consider another example

Solve the equation for t from 0 to 8 :

dx/dt =   x -  x y
dy/dt = - y +  x y

Initial condition xt=0 = 0 , yt=0 = 2

To rename the variables, t is set to be 't0'. As the first differential equation is dx/dt, we set x to be 'x0'. The second differential equation is dy/dt, we set y as 'x1'. The renamed problem is :

Solve the equation for t from 0 to 8 :

dx0/dt0 =   x0 -  x0 x1
dx1/dt0 = - x1 +  x0 x1

Initial condition x0t0=0 = 0 , x1t0=0 = 2

Input 4 : a='x0-x0*x1'
 

Input 5 : b='-x1+x0*x1'
 

Input 6 : r=rkode('a|b','1|2',0,8,0.1)
 

The resulted r is a matrix consists of the elements. The first is t0. The second is a matrix consists of two lists. One list is the solution for x0 and the other is for x1.
To plot both x0 and x1 against t0 :

Input 7 : d=r(1)
 

Input 8 : plot('line','r(0):d(0)|r(0):d(1)',0,8,0.1)
 

The plot is shown below :

 To see more clearly how the two variables x0 and x1 are related, we plot a graph of x0 against x1 :

Input 9 : plot('line','d(0):d(1)',0,8,0.1)
 

The values of x0 and x1 at any time are at the perimeter of a circle-liked loop. This example is actually a model for 'predator-prey' in ecology and x0 and x1 are the populations of the predator and the prey and they follow the a closed loop of change, meaning that the populations of the two species are changed in synchronization and the changes is repeated periodically.
Ordinary differential equation with second order derivative can also be solved by 'rkode'. Below is an example demonstrates this :

Solve the equation for t from 0 to 10 :

d2x/d2t + 10x = step(t-3) where 

step(t-3) = 1 if t >= 3
          = 0 if t <  3 

Initial condition xt=0 = 0, dx/dtt=0 = 1

Here we rename t as 't0'. x as 'x0'. But in this case, we also set a new variable 'x1' to be dx/dt. With this change, the original second order differential equation would separate into two first order differential equations :

Solve the equation for t from 0 to 10 :

dx0/dt0 = x1
dx1/dt0 = step(t0-3) - 10 x0

where,

step(t-3) = 1 if t >= 3
          = 0 if t <  3

Initial conditions : x0t0=0 = 0, x1t0=0 = 1

Input 10 : r=rkode('x1|step(t0-3)-10*x0','0|1',0,10,0.1)
 

To plot x0 against t0,

Input 11 : d=r(1)
 

Input 12 : plot('line','r(0):d(0)',0,10,0.1)
 

The plot is shown below :

The above example is actually the equation describing a damped oscillation of a particle attached to a spring with an external constant force equal to one acting on it when t >= 3.