HMath-4010: Solving Ordinary Differential Equations

The first part of this tutorial considers solving ordinary differential equations (ODE) in HyperMath. A second order differential equation is used for illustration purposes as they are more common. The equation can represent a mechanical spring-mass-damper system or a series RLC circuit, both excited by an arbitrary sinusoidal input of amplitude A.

For the spring-mass-damper system:

Φ (referred to as the state variable) represents position.

α damping coefficient divided my mass.

β the spring stiffness divided by the mass.

The excitation is a force.

For a series RLC circuit system:

Ф represents current.

α resistance divided by the inductance.

β the elastance divided by the inductance.

The excitation is a voltage.

Step 1: Converting to first order equations.

1.	This is a required step since the ODE solvers can only solve first order explicit equations. This is done by introducing one additional state variable, ξ, to be equal to the derivative of the original state variable, Ф. Now, there are two coupled first order differential equations:

The new state variable relationship equation:

And the original equation, in terms of the new state variable:

Note that the equations are expressed in explicit form (only derivatives on the left side). The task is to solve these two coupled equations for the two state variables as a function of time.

To simplify, we will assume that all the parameters (α, β, A, ω) are equal to 1 without any loss of generality.

Step 2: Implementing the two equations.

1.	The two equations derived in the above step need to be implemented inside a function. For this tutorial, we will name the function ODE_EQN. The function needs to have the following characteristics:

•	Take the current time step point as the first input argument. It is a scalar.

•	Take the current values of the state variables as the second input argument. It is a column vector of length matching the number of state variables in the system (two in our case).

•	It must output the derivatives of the state variables in a column vector (in our case, the two equations). The order of the derivatives in the output is arbitrary, but must match the initial conditions when the ODE solver is called.

This is what our function looks like:

function ODE_EQN (t, x)

{

// t is a scalar

// x is a column vector of 2 elements, representing the states Ф and ξ respectively

// set all parameters to 1

A = 1; alpha = 1; beta = 1; omega = 1;

// implement the equations next

dPHI_dt = x(2); // first equation

dXI_dt = A*Sin(omega * t) - alpha*x(2) - beta*x(1); // second equation

out = [dPHI_dt; dXI_dt]; // output as a column vector

return out;

}

Step 3: Using the ODE solver.

1.	With the equations defined in our own function, the next task is to use a HyperMath built-in ODE solver. For this tutorial, we will use the RK45 (Runge-Kutta mixed fourth-fifth order) function, specifically in its following form:

x = RK45(func, init, time, reltol, abstol, userdata)

where

Option	Description
func	The name of the function we just created.
init	The initial conditions of the state variables, specified as a vector.
time	The time point where the values of the two state variables are to be calculated.
reltol	Relative tolerance of convergence. We will use the default value of 1e-03 and omit this.
abstol	Absolute tolerance below which the solutions are unimportant. We will also use the default 1e-06.
userdata	This is not required for our example and will be omitted.
x	A matrix of the solutions at each time point. The columns would represent the state variables.

The initial conditions and the time points where we would like to evaluate the solutions need to be defined. For this example, we would assume the states start at zero values and create a row vector of two zeros to represent these (the order of state variable representation in this vector must match the order assumed in the second argument of the function we created). We choose to evaluate the solutions starting from time zero to ten in increments of 0.1.

With the above two defined, we can simply call RK45 with our function name, initial condition, and time points. Then, the solver starts the numerical integration starting with the initial values. It calls the function at each time step, passing in the current time and state variable values to obtain the corresponding derivatives from our function to calculate the next time step values until it reaches the last time point. This all happens behind the scene.

The above can be implemented as follows:

initial = [0, 0];

// Initial conditions for the states

t = [0:0.1:10];

// The time steps for solution

x = RK45("ODE_EQN", initial, t);

The output variable x contains the time-series values for both states. It has two columns and 101 rows. Each column represents a state variable, the first one being Ф and the other ξ (consistent with the order assumed in the user function). Each row represents a time-step.

We can extract them and plot against time as follows:

y1 = x(:,1); // Extract the first state response

y2 = x(:,2); // Extract the second state response

PlotLine(t, y1,'new');

SetLegend("Phi");

PlotLine(t, y2);

SetLegend("Xi");

SetTitle("RK45 Solution");

For spring-mass-damper system, the plot represents the position and velocity. For the RLC circuit, it represents the current and its rate of change.

HMath-4010: Solving Ordinary Differential Equations

HMath-4010: Solving Ordinary Differential Equations

HMath-4010: Solving Ordinary Differential Equations

See Also: