A comparative analysis of different approaches used for modelling and solving differential equations in Simulink/MATLAB

2.1. Solving differential equations with non-zero initial conditions

A differential equation is a mathematical equation for an unknown function of one or several variables that relates the value of the function itself and its derivatives of various orders. Ordinary differential equations (ODE) are derivatives of a dependent variable with respect to one independent variable, usually time. Ordinary differential equations (ODE) play a fundamental role in the modern technological systems where circuit and systems can be described in different ways using differential equations of first order (e.g. RC and RL), second order and higher order for analyzing their dynamic performance. Second-order systems (mass-spring-damper systems, RLC circuits, etc.) are encountered frequently in practical applications [10] and needed to be solved for analysis. To simplify analysis, higher-order systems are approximated as second-order systems.

Thus, a lot of electric systems are modeled mathematically using different approaches such as a differential equation, or by a transfer function and state space for determining a change in the behavior of a system over time [11-13]. These approaches are discussed below.

2.2. Integrator based approach

It includes defining systems (related to fields such as communications, controls, signal processing, video processing, image processing, etc.) by linear equations, or differential equations, or difference equations. Afterwards, implementing these equations with graphical block diagrams that include integrator blocks available in Simulink library.

Simulink allows simulation of linear time-invariant (even nonlinear or time-varying) differential equations that can be written in an explicit form given by:

The corresponding block-scheme (to be designed in Simulink) is as shown in Fig. 1.

Fig. 1The integrator block scheme for solving differential equations in Simulink

In order to determine a linear dynamic system’s time response that depends on initial conditions and the system input, a set of mathematical equations representing the system is derived. For some simple systems, a closed-form analytical solution can be easily obtained. But for most of the systems (including nonlinear systems or those subject to complicated inputs), an integration must be carried out numerically to determine the output response that represents how the state of a dynamic system changes in time when subjected to a particular input. Simulink/MATLAB provides all useful resources for calculating time responses to several types of inputs.

As the derived equations consist of differential equations, an integration is required to be performed in order to determine a time response of the system. The integrator block with 1/$s$ as the block symbol available in Simulink library is used to calculate the integral of the signal provided at the input i.e. left side of the block [14-15]. While solving a differential equation, the major task to perform includes integrating a given equation so as to get back the function without employing the derivative function. Integration of the derivative of a function is equal to the function itself. The integrator block serves the same purpose in the Simulink model. Simulink integrator blocks used in a block diagram are equal to an order of the differential equation that is being analyzed i.e. for solving a first order differential equation, only one integral block is needed, whereas for a second order differential equation, the number of blocks dragged and dropped into the model is two.

2.3. Transfer function approach

Simulink can also solve ordinary differential equations by defining systems using transfer functions, and thereby using techniques for numerical integration. Simulink has a standard block called the Transfer Fcn block which allows the direct input of numerator and denominator terms of the transfer function describing a given dynamic system. The Transfer Fcn block has a default form of first-order in the denominator ($s+1$) which can be modified for the specified different order for the numerator and denominator polynomials by entering coefficients associated with the respective polynomials (zero coefficient of the missing term must be included).

Any linear system can be characterized using a transfer function $G$, which is an operator that takes an input $u$ to the output $y$. For a continuous time system, system transfer function $G$ is essentially a ratio between Laplace transform of the output response of a system $Y\left(s\right)$ and Laplace transform of the input forcing function $U\left(s\right)$ i.e. $G\left(s\right)=Y\left(S\right)/U\left(S\right)$. The transfer function for any system can be obtained by taking a Laplace transform of the differential equation representing the given system with all initial conditions set to zero, and then solving the output response in terms of a given input forcing function $u\left(t\right)$ with a factor multiplying $L\left[u\left(t\right)\right]$ taken out as the desired transfer function [2]. The desired transfer function is defined in Simulink by including coefficients of the polynomials in numerator and denominator of the transfer function as the inputs to the Transfer Fcn block in the Num and Den fields, respectively. This can be done by accessing the Function Block Parameter box where all information about the system can be entered within one block i.e. the Transfer Fcn block.

Now, the target is how to solve differential equations having non-zero initial conditions using transfer function approach. The paper presents two ways for resolving this issue (Section 3).

2.4. State space approach

A system in its state space representation replaces an $n$th order differential equation with a single first-order matrix differential equation. The method is convenient for breaching a higher-order differential equation into sequences of first-order equations so that they can be easily solved using a matrix method. In order to model a second-order ODE using state space, the equation is broken down into two first-order state equations [16-18]. This can be achieved by allocating a subscripted variable for each state of the system in an order of increasing derivatives which includes letting $x_1=x$, $x_2=x^{\text{'}}=x_1^{\text{'}}$, and$x_2^{\text{'}}=x\text{'}\text{'}$. This allows a re-statement of the equation in terms of state variables only.

To begin, select the State-Space block from the Continuous sub-menu of the Simulink library.

At this point the model is very general, and an equation of any order can be set up for solution in the block parameters. The equation inside the State-Space block is $x^{\text{'}}=Ax+Bu$, $y=Cx+Du$. This represents the basic state-space equation, where $x^{\text{'}}$ = a vector of the first-order state variables, $y$ = output vector, $x$ = the state variable vector, $u$ = a forcing input function, $A$ = the state matrix, $B$ = input matrix, $C$ = output matrix, and $D$ = transmission matrix.

3.1. Integrator method

For solving Eq. (4) using integrator based approach, first step is to rearrange the given differential equation. Keep the highest order derivative on the left-hand side, and move all the remaining terms on the right-hand side. The resulting equation can be simulated in Simulink as shown in Fig. 3. The initial conditions specified at zero time instant i.e. $t=$ 0 can be entered easily in the Integrator block.

Fig. 3The Simulink circuit model using integrator and transfer function (method 1)

Simulink (a MATLAB component) is very useful for simulating dynamic systems using a block-diagram model based approach [19]. A block diagram is simply a graphical representation of a process (which is composed of an input, the system, and an output). Thus, an interconnection of blocks in Simulink represents the system’s mathematical model. The lines interconnecting the blocks represent the variables describing the system behavior. The Display Sink block is a digital readout of an input signal at a current simulation time (the specified simulation stop time). The Scope Sink Block is used to display a signal as a function of time in a Simulink model. The term model refers to a set of mathematical equations representing a physical system, and thereby relates the system’s output signal to its input signal.

The simulated output after running simulation from $t=$ 0 to $t=$ 10 is shown below in Fig. 4. The Display Sink block is showing the numerical output value of IVP given by Eq. (4) at a time instant equal to 10 s (the specified current simulation stop time). Sampling issues must be taken care of while running a simulation in Simulink (a front-end tool designed for solving ordinary differential equations numerically). If the variable-step solver option is selected, the maximum step size parameter in the Max Step Size field box present in the Configuration Parameters window should be set to a lower value for better simulation results.

Fig. 4The simulated output using integrator method

3.2. Transfer function method

For using transfer function approach, take Laplace of Eq. (4). While writing the Laplace, following properties should be used:

The Laplace transform of Eq. (4) after plugging in the initial conditions becomes:

After solving for $Y\left(s\right)$ and combining into a single term, gives:

or

The analytical solution can be computed by applying partial fraction method on Eq. (8), and then taking inverse Laplace transform of the terms obtained after doing partial fraction decomposition [20]. The exact solution comes out to be:

Now, the target is how to solve it digitally using circuit simulation.

3.3. Method 1

3.3.1. The newly developed first way for solving Eq. (8) using transfer function in Simulink (method 1)

Take the factor multiplying the Laplace of an input function $2\sin \left(3t\right)$ i.e. the multiple of $6/s^2+9$ as the transfer function. Thus, the transfer function to be entered in a Simulink Transfer Fcn block becomes equal to as follows:

11

$Y\left(s\right)=\frac{\left({-s}^3+{2s}^2-9s+24\right)}{\left(6s^2-36s+90\right)}L(2\sin \left(3t\right).$

The above transfer function can be re-written as:

12

$Y\left(s\right)=\frac{s\left({-s}^3+{2s}^2-9s+24\right)}{\left(6s^3-36s^2+90s\right)}L(2\sin \left(3t\right).$

The transfer function $s\left({-s}^3+{2s}^2-9s+24\right)/\left(6s^3-36s^2+90s\right)$ can be implemented in Simulink by using a Derivative block in series with a Transfer Fcn block as shown in Fig. 3. The simulated output result of the above stated transfer function (method 1) by implementing Eq. (12) is shown below in Fig. 5.

Fig. 5The Simulated output with the transfer Fcn block using the method 1

3.4. Method 2

3.4.1. The newly developed second way for solving Eq. (8) in Simulink (method 2)

Take the factor multiplying the Laplace of a unit step function$i.e.\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{}\mathrm{o}\mathrm{f}\mathrm{}1/\mathrm{s}$ as the transfer function. Now, implement the transfer function $\frac{\left({-s}^4+{2s}^3-9s^2+24s\right)}{\left(s^4-6s^3+24s^2-54s+135\right)}$ as shown in Fig. 8.

Fig. 6The simulated output values for time instants t= 0 s and t= 10 s using an integrator scheme

Fig. 7Outputs for t= 0 s and t= 10 s using the Method 1 (based on transfer function)

Fig. 8The Simulated model with the transfer Fcn block employing the method 2

By default, the Step block’s step time (time, in seconds, when an output steps/changes from the Initial value parameter to the Final value parameter) is 1. Also, the Initial value parameter (the block output until the simulation time reaches the Step time parameter) is 0 by default whereas, the Final value parameter (the block output when the simulation time reaches and exceeds the Step time parameter) is 1 by default.

Fig. 9 shows the output using a second way of employing transfer function approach (method 2).

The simulated outputs for time instants $t=$ 0 s and $t=$ 10 s using integrator approach are –1 and –7.366×10¹², respectively.

Fig. 9The simulated output with the transfer Fcn block using the method 2

3.4.2. Exact output

The analytical solution given by Eq. (10) can be simulated as shown in Fig. 10.

Fig. 10The simulated model for generating an exact output of Eq. (10) using analytical solution

The Clock Source Block is used in the Simulink model shown in Fig. 10 for generating a signal equal to the current time in the simulation. The Simulated exact output describing Eq. (10) is shown in Fig. 11.

Fig. 11The simulated exact output using an analytical solution

3.4.3. Output using state variable approach

Select the state variables for the Eq. (1). Choosing the state variables as successive derivatives, we get two first order differential equations as shown below:

13

$x_1^{\text{'}}=x_2,$

14

$x_2^{\text{'}}=-15x_1+6x_2+2\sin \left(3t\right)=-15x_1+6x_2+f\left(t\right).$

These can be reorganized to obtain:

15

$\left\{\begin{array}{c}{x}_1^{\text{'}}\\ x_2^{\text{'}}\end{array}\right\}=\left[\begin{array}{cc}0& 1\\ -15& 6\end{array}\right]\left\{\begin{array}{c}{x}_1\\ x_2\end{array}\right\}+\left[\begin{array}{c}0\\ 1\end{array}\right]f\left(t\right),$

16

$y=x=x_1=\left[\begin{array}{cc}1& 0\end{array}\right]\left\{\begin{array}{c}{x}_1\\ x_2\end{array}\right\}+\left[0\right]f\left(t\right).$

The Simulink model for simulating Eq. (15-16) using the State Space block is shown in Fig. 12. The block allows entering the direct input values of the A, B, C and D matrixes that are unique to a particular state-space model by going into the block parameter dialog box. The initial conditions are specified separately in a place to enter the ICs in the State Space block parameter window.

Fig. 12The model for comparing the state space solution with the exact solution

The model in state-space can be easily mapped in Simulink environment for numerical simulation. The simulated output using state space approach is presented in Fig. 13.

A comparison of the exact output with the state space based output is shown in Fig. 14. A spirit of comparing the simulation result with an analytical one is to ensure that a numerical approximation gives an acceptable result.

Fig. 13The simulated output using state space approach

Fig. 14Comparison of the simulated exact output with the state space based output

The Simulink model for comparing the outputs using an analytical approach, state space method, integrator method and transfer function approach (method 2) is shown in Fig. 15.

Fig. 15The model for comparing the outputs based on different approaches

A comparison of outputs using an analytical method, state space method, integrator method and transfer function (method 2) using the analytical solution as a benchmark is shown in Fig. 16.

Fig. 16Comparing the simulated outputs