Comparative Study of Matlab ODE Solvers for the Korakianitis and Shi Model

Changing parameters of the Korakianitis and Shi heart valve model over a cardiac cycle has led to the investigation of appropriate numerical technique(s) for good speed and accuracy. Two sets of parameters were selected for the numerical test. For the seven MATLAB ODE solvers, the computed results, computational cost and execution time were observed for varied error tolerance and initial time steps. The results were evaluated with descriptive statistics; the Pearson correlation and ANOVA at α0.05. The dependence of the computed result, accuracy of the method, computational cost and execution time of all the solvers, on relative tolerance and initial time steps were ascertained. Our findings provide important information that can be useful for selecting a MATLAB ODE solver suitable for differential equation with time varying parameters and changing stiffness properties.

and Correct, Evaluate (PECE) method [21]. Adams' formulas have successfully been used to solve many differential equations [32]. Spijker investigated their stability and convergence property [2] while their effectiveness for treating stiff problems was reported by Lavernano and Mazzia [33].
The ode23t solver implements the trapezoidal rule. The trapezoidal rule is a one-step low order but absolutely stable technique [26]. It is the simplest implicit technique for solving ODE [34] but has a low efficiency [22]. The ode23s solver is a more efficient one step solver in MATLAB ODE suite. It implements the modified Rosenbrock techniques. The general form of Rosenbrock techniques is discussed in the literature [35]. The modified versions [9,36] have improved stability property. A very reliable solver for stiff problems is the Backward Differential Formulas (BDFs). Their development and implementation are discussed [2]. The numerical differential formulas (NDFs) are a modification of the BDFs [23]. Being modified versions of BDFs, NDFs give better results but cannot be implemented beyond order six as they are unstable at higher order [37]. The MATLAB ode15s solver implements the NDF and BDF alternatively, while the ode23tb implements the trapezoidal rule as a first stage and then the BDF formula as a second.
Omale et al. reported the solution to six sets of standards equation with the MATLAB ODE solvers [38]. Abelman and Patidar [39] compared the numerical results of three MATLAB ODE solvers (ode45, ode15s and ode23s) with that of the nonstandard finite difference method (NSFDM). They showed that the NSFDM can be used to obtain good result without step size restrictions. Yatim et al [40] compared results they obtained using a newly developed backward Differential formula with those obtained using ode15s and ode23s for three sets of problems. Ocherova [41] modified the ode23 solver and showed that the relative tolerance has effect on the solution, evaluation time, computation cost and accuracy of the solver. Dallas and Machairas [20] investigated the MATLAB solvers and their parameters for a class of nonlinear hybrid system. He attempted to determine which of the solver was best in terms of speed and accuracy for such type of problems.
The aim of this paper is to evaluate the effect of tolerated error and initial time step on the accuracy, execution time and computational cost of the seven MATLAB ODE solvers for an equation with time varying parameters that may become stiff at certain parameter levels. The other sections of this paper are organized as follows: Section two contains a brief description of the KS model and the selected parameters and variables for simulation. Section three of the paper presents the analysis and discussion of results. The paper ends with conclusion and remarks based on the findings.

The KS Heart Valve Model
The KS heart valve model [3] that describes the dynamics of the heart valve is shown in equation (1) below: represents the angular displacement of a heart valve. 1 is the pressure of the chamber or vessel from which blood is flowing through the valve while 2 is the heart chamber or blood vessel which blood flows into from the valve. Consequently, ( 1 − 2 ) is the pressure gradient across the valve. Q is the flow rate through the valve. K p , K f , K b and K v are time invariant parameters representing the pressure, friction, blood flow and vortex constant respectively [3]. Their values for the aortic and mitral valve respectively are given in Table 1: 3.5 rad/(s m) Source: [3] Equation (2) is a simplified form of equation (1) above This can be represented with a system of two 1 st order ordinary differential equations as follows:

Model Parameters for Simulation
Prior simulation was used to observe the dynamical behavior of the model on a time series plot. After then, two sets of parameter values of A, B and C that yielded differing stability property were selected for this study. They are presented as follows:

Numerical Experiments
For the numerical experiments, the valve was assumed to be opening from its rest position, so initial values of (0, 0) was used for ( 1 and 2). An average time over which the semilunar valves open is 0.3 Seconds [42][43]. Thus, selected iteration time for the simulation is 0.3 Seconds. A MATLAB function file for equations 2a and 2b was written for Simulations purpose. While varying the relative tolerance and initial time steps of the schemes, the following were observed for the seven (7) MATLAB ODE solvers: 1. The calculated results for angular displacements.  Fig. 1 shows the time series plot of angular displacement for the KS model over an iteration time of 0.3 Seconds for both cases of parameter set. As observed, the solution is a decaying sinusoid that settles to about 89.3978degrees for case 1 while the solution decays exponentially and settles to about 73.3152 degrees for case 2. The change in dynamics of the system is due to the change of stability properties that occurs as parameter changes. As observed in Table 2, the calculated values of angular displacement for a very low relative tolerance was unaffected by the type of MATLAB solver as the value of P=0.3846 is greater than 0.05. The values of the angular displacements for the two different sets of parameter (Case 1 and Case 2) vary significantly as P=1.71E-15 is far less than 0.05.  The average angular displacements, standard deviation and their correlation with relative tolerance, for the different solvers, are shown in Table 3. For all the solvers, there is a correlation between the calculated angular displacements and the relative tolerance. However, the correlation coefficient changes from negative to positive or vice versa for the two cases of selected parameters, for all solvers except ode113. The values of angular displacement for ode45 and ode23tb solvers were least affected by the tolerated errors, as they have the least standard deviations. Table 4 is the ANOVA table for case 1 parameters. Similar results were obtained for case 2 parameters. As shown, the P-values for rows and columns were 0.4592 and 0.4068 respectively, both being greater than 0.05. This is an indication that both the type of solver and the relative tolerance do not significantly affect the result of angular displacements obtained for the selected initial time step.

34
Volume 19 The average angular displacements, standard deviation and their correlation with initial time step, of the different solvers, are shown in Table 5. For all the solvers, there was a strong correlation between the calculated angular displacements and the initial time step. The correlation of ode23t does not vary much with the model parameters. All other schemes, but ode113, had correlation coefficients which changed from negative to positive or vice versa for the two cases of selected parameters. The deviation of angular displacements for all solvers are very small indicating that the calculated result of a solver was not very much affected by the initial time step. Table 6 shows the ANOVA table for case 1 parameters. Similar results were obtained for case 2 parameters. As shown, the P-values for rows and columns were 3.24E-42 and 0.4759 respectively. The calculated angular displacement was affected by the type of solvers for the selected relative tolerance, while the result of angular displacement is unaffected by the initial time step (P<<0.05 for columns). The number of times the integrand is evaluated defines the computational cost [37]. As expected, computation costs were higher for parameters of case 1 than that of case 2 due to higher degree of non-linearity. As observed in Figure 2, the computational costs are higher for low relative tolerance for all the solvers. Table 7 shows the ANOVA table for case 1 parameters. Similar results were obtained for case 2 parameters. The P-values for rows and columns were 0.002213 and 0.000535 respectively, indicating that the computational cost was affected by the type of solvers and the relative tolerance selected. (P<<0.05 for both rows and columns).  Fig . 3 shows the log values of computational costs for varied initial time steps. The ode23s is the most expensive technique, followed by the ode23tb and ode23s respectively. The implicit solvers are more expensive than the explicit ones because of the evaluation of Jacobian matrices and derivation of partial derivatives in its procedure for improved efficiency [23]. The multistep techniques are a function of already computed values of preceding steps and may need to evaluate the differential function a minimum number of times. Consequently, ode113 and 0de15s are the least expensive. Although ode45 and ode23 are both Runge-Kutta scheme, ode45 is less expensive because it is a higher order Runge-Kutta scheme over the ode23. Table 8 shows the ANOVA table for case 1 parameters. Similar results were obtained for case 2 parameters. The P-values for rows and columns were 1.83E-07 and 0.527274 respectively, indicating that the computational cost was affected by the type of solvers (P<<0.05) but unaffected by the initial time step (P>0.05 for columns).   4 shows the execution time for the different solvers for different relative tolerance. The execution time is expected to be related to the computation cost of a solver. So as expected, the implicit schemes had higher execution time, with ode23s being the highest followed by ode23tb and ode23t respectively. The explicit multi-steps technique (ode15s) had the least execution time, followed by the Runge-Kutta techniques (ode45 and ode23).  Table 9 shows the ANOVA table for case 1 parameters. Similar results were obtained for case 2 parameters. The P-values for rows and columns were 0.001548 and 0.00038 respectively, indicating that both the solver and relative tolerance affected the execution time (P>0.05 for columns).    Table 8 shows the ANOVA table for case 1 parameters. Similar results were obtained for case 2 parameters. The P-values for rows and columns were 9.5E-29 and 0.025105 respectively, indicating that the execution time was affected by the type of solvers and the preset initial time steps (P<0.05 for both cases).  The percentage number of unsuccessful steps in Table 11 shows the percentage of times a step was repeated due to re-computation of the initial time step. It is related to the local error generated by the respective method embedded in the solver. As observed, the percentage failed steps has positive correlation with the selected relative tolerance. It implies that the smaller the relative tolerance, the lower the local computational error. The ode23tb is the most accurate solver regardless of relative tolerance and parameter combination, having the least number of failed steps.   Fig. 6 shows the percentage of unsuccessful steps of the different solvers for varied initial time steps. The percentage of failed number of steps is highest for the ode15s scheme, followed by the ode45 and ode113 schemes regardless of the initial time steps and parameter set. The implicit solvers had the least number of failed steps, indicating that they are more accurate schemes. Table 13 is the ANOVA table for case 1 parameters. Similar results were obtained for case 2 parameters. The P-values for rows and columns were 4.85E-14 and 4.22E-06 respectively, indicating that the execution time depends on both the solver and the initial time step (P>0.05 for columns).

Discussions and Concluding Remarks
The aim of this work was to evaluate the performance of the seven MATLAB ODE solvers for the Korakianitis and Shi model with changing parameters. The choice of error tolerance is fundamental to the accuracy and cost of results regardless of the type of solver, whereas, the effect of the initial time step is insignificant. However, the value of initial time step can affect the number of unsuccessful steps, and hence the execution time. If computation time is a priority, the multistep techniques (ode113 and ode15s) are best suitable as they are the least expensive. If a high level of accuracy is desired, the implicit solvers, ode23tb followed by ode23s and ode23t respectively are the best alternatives. However, the high accuracy of the implicit solvers will be a tradeoff for computation time as they are relatively costly. For a non-stiff equation, the ode45 solver is a good choice for first alternative as it presents a balance between accuracy and computational cost.