# Mathematical Modeling of Real-Time Systems Using Heun and Piecewise Methods

Authors: Urfa Malik Gul; Anand Paul; K.-W.-A. Chee
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

It is often said that mathematical modeling is an implementation of mathematics in real-world problems with the aim of better understanding them, so we can say that mathematical modeling is linked to the solution of problems. Some of the essential principles and procedures of mathematical modeling are discussed using formulas and equations. We investigate the stability and convergence characteristics and demonstrate the suitability of different mathematical methods in a set of numerical examples. The described methods in our paper are the best choices for the simulation of linear phenomena and are more efficient for use with high-order spatial discretization. We emphasized the importance of mathematical modeling technologies used in computational tools. Our study shows that these new methods are more stable with lower errors.


---

## Body

## 1. Introduction

Various numerical methods are used to obtain the initial value solution, and studies show great development on this topic in this century. Other factors also amount to the reason for the implementation of mathematical methods like electronic computers needing efficient numerical algorithms in mathematical modeling. Many complex methods that are driven by partial differential equations (PDE) evolve [1–3], which requires the search for the appropriate numerical solution to a system of ordinary differential equations (ODEs) in a wide range of situations [4–6]. The differential equations contain derivatives, either ordinary or partial derivatives, and along with an initial condition that determines the value of the unknown function at a certain point, make up the initial value problem (IVP). These equations are also used in semiconductor engineering, biology, chemistry, physics, and economics to solve various problems in research. There are several analytical approaches for determining the solution of differential equations [7, 8]. On the other hand, analytical methods cannot always solve complex differential equations. Therefore, it is important to make numerical approximations.The advanced differential equation is solved by using numerical techniques for which computer programming is a very powerful tool. The well-known Runge–Kutta (RK) analysis is the most often used technique for integrating these technologies of ODEs in real time [7–10]. RK methods can be classified in a variety of ways, such as the structure of implicit or explicit, by convergence order, or the number of steps [1]. The complexity of the equation that must be solved determines the numerical stability of various numerical methods.Figure1 summarizes the flow of ordinary equations used to model real-time systems and the application of analytical and numerical methods to solve these equations. An example here is encapsulating the physical characteristics of a parachutist descending, where g is gravitational constant, m is mass, c is drag coefficient, and dv/dt is an unknown function/differential equation (rate equation).Figure 1
The summary of steps in the application of ODEs to engineering problem solving, such as a parachutist descending.In the RK family, the most well-known member is the classic Runge–Kutta method (RK4), and (1) presents the initial problem as follows:(1)dydt=ft,y,yt0=y0.In the above equation (1), y can be a scalar or vector function that is not known at time t that we want to estimate. dy/dt is the time rate at which y changes and t0 is the initial time; hence, the corresponding y value is y0.(2)yn+1=yn+12hk1+2k2,tn+1=tn+h,k1=ftn,yn,k2=ftn+h2,yn+hk12,k2=ftn+h,yn+hk1.Thek1 and k2 are known as stages of the RK method which are shown in (2) and (3). They relate to various estimations of solutions for the slope. We note here that yn+hk1 is the Euler step from tn,yn with stepsize h. Now, the following equation shows an explicit second-order RK method by supposing(3)yt+h=yt+hb1k˜+b2k˜2+Oh3,k˜1=ft,y,k˜2=ft+c2h,y+ha21k˜1.Therefore, the second-order RK supposition becomes(4)ft+c2h,y+ha21k˜1=ft,yc2hftt,y+ha21fyt,yk˜1+Oh2,ft,yc2hftt,y+ha21fyt,yk˜1ft,y+Oh2,yt+h=yt+hb1ft,y+b2ftt,y+c2hftt,yha21fyt,yft,y+Oh3,yt+b1+b2hft,y+b2h2c2hftt,y+a21fyt,yft,y+Oh3.Therefore, for four unknowns, we have three nonlinear equations which lead to Euler’s method (first-order method). Euler’s method is commonly used to build more complicated approaches like the predictor/corrector method and is sometimes also called the midpoint rule as follows:(5)yn+1=yn+hk2,k1=ftn,yn,k2=ftn+h2,yn+h2k1.As we know, the RK methods are a set of implicit and explicit optimization methods in numerical analysis, including the well-known Euler’s method, which is usually used in temporal discretization for approximate solutions of ODEs [1, 11]. Here, the differential equation is considered as a formula that can be used to compute the slope of the tangent line to the curve at any point once we determined the location of that point. The solution to the ODE at a time tn:yn≈ytn is approximated by the value of yn. The explicit Euler’s method shows that the solution yn+1 is an explicit function of yi for i≤n. We must mention here that there is no assurance about any curves’ concavity that will remain consistent. It may shift from concave-up to concave-down. Because the real solution curve goes up and down, there is no approximate way to estimate whether our next point will be an over- or underestimate. To overcome this deficiency, we have a clever Heun’s method to use.


## 2. Heun’s Method

Heun’s method is the improved (or modified) Euler method and the second order of the RK method. The modified Euler method allows us to find a clear expression fory given a finite number of elementary functions, x. It is possible to solve the ODE for y and x given their initial values.The Euler’s method algorithm is as follows:(1)
Start(2)
Define function(3)
Determinex0, y0, h, and xn (here x0 and y0 are initial conditions while h is interval and xn is required value)(4)
n = (xn − x0/h+ 1)(5)
Fromi = 1 ton start looping(6)
y = y0 + h × f (x0, y0), x = x + h(7)
Print they0 and x0 values(8)
Ifx < xn, assign x0 = x and y0 = yIf no, go to 9(9)
Loopi end(10)
StopWe can use the algorithm and flowchart in a high-level programming language to construct a program for Euler’s method though this approach is not one of the best techniques to solve ODEs.Heun’s formula is used to solve the ODE which is given in the initial state [10, 12]. Euler’s method is the foundation of Heun’s method. As we discussed earlier, Euler’s method fails to converge in the case of a concave-up curve because it overestimates the next point. Heun’s method solves this solution problem by taking two tangents on each side of the curve. One tangent underestimates and the other overestimates, and then, Heun’s method estimates the next point from both tangents using Euler’s method. Using (6), we can determine the numerical solution to the problem of initial values using the following method:(6)yi+1=yi+k1+k2h2,k1=fxi,yi,k2=fxi+h,yi+k1h.where k1andk2 are under- and overestimates. The step size is denoted by h, and initial conditions are xiandyi. The accuracy of Euler’s method improves linearly as step size decreases but the accuracy of Heun’s method improves quadratically [10, 12–15].


## 3. Derivations

While using the elementary idea that the line’s slope equals its rise over run, the parameters at the endpoint can be calculated using the following symbolic development:(7)slopeleft=fxi,yi,sloperight=fxi+h,yi+hfxi,yi,slopeideal=Δyh,slopeideal=12slopeleft+sloperight,Δy=hslopeideal,xi+1=xi+h,yi+1=yi+Δy,yi+1=yi+hslopeideal,yi+1=yi+12hslopeleft+sloperight,yi+1=yi+h2fxi,yi+fxi+h,yi+hfxi,yi.Using Euler’s method, the next point in the numerical solution is estimated roughly, and the initial estimation could be predicted or corrected using this information. Therefore, we can calculate the slope of the prediction line by taking the average of the slopes of the left and right tangent lines, at either end of the interval. The values of thef (x, y) on the right-hand side of Equation (7) can be determined [16]. Following a rough estimate of the location of the next solution point from Euler’s method, the coordinates may be used to determine the slope of the tangent line at the right end of the interval.


### 3.1. Modified Algorithm

The initial or boundary conditions must be given for both the algorithm and the flowchart, as follows:(1)
Start(2)
Define function for slope calculationf (x, y)(3)
State variables(4)
Input value(5)
Find slope while using initial value(6)
Find newy, y = y0 + m0 × h(7)
Increase value ofx, which is x1 = x0 + h(8)
Find new slope usingx1 and y(9)
Calculate the mean ofm0 and m10)
(Findy new and assign y new with y again(11)
From step 6, repeat until two consecutivey are equal(12)
Repeat step 5 untilx = xn(13)
Printx and y corresponding(14)
Stop


### 3.2. Flowchart

The numerical solution of ODEs with Heun’s method is simple from both the numerical and programming points of view as compared to analytical techniques. The updated Euler’s method algorithm and flowchart are described in Figure2. They may be used to produce source code in any high-level programming language.Figure 2
Modified Euler’s method flowchart.In the case of complex problems, analytical methods frequently fail; however, Heun’s method does not fail but provides better solution accuracy and high speed. It also improves Euler’s technique with a lengthy code where we can pick a lower number ofh to get high accuracy.


## 3.1. Modified Algorithm

The initial or boundary conditions must be given for both the algorithm and the flowchart, as follows:(1)
Start(2)
Define function for slope calculationf (x, y)(3)
State variables(4)
Input value(5)
Find slope while using initial value(6)
Find newy, y = y0 + m0 × h(7)
Increase value ofx, which is x1 = x0 + h(8)
Find new slope usingx1 and y(9)
Calculate the mean ofm0 and m10)
(Findy new and assign y new with y again(11)
From step 6, repeat until two consecutivey are equal(12)
Repeat step 5 untilx = xn(13)
Printx and y corresponding(14)
Stop


## 3.2. Flowchart

The numerical solution of ODEs with Heun’s method is simple from both the numerical and programming points of view as compared to analytical techniques. The updated Euler’s method algorithm and flowchart are described in Figure2. They may be used to produce source code in any high-level programming language.Figure 2
Modified Euler’s method flowchart.In the case of complex problems, analytical methods frequently fail; however, Heun’s method does not fail but provides better solution accuracy and high speed. It also improves Euler’s technique with a lengthy code where we can pick a lower number ofh to get high accuracy.


## 4. Piecewise Method

Nonlinear programming is used to solve a variety of optimization issues in science and technology. Piecewise linearization methods are used in recent years to turn nonlinear programming into linear programming or a mixed convex programming problem to get the desired solution [16–19]. The piecewise function is a hybrid function that is defined by numerous subfunctions and each corresponding to a different interval of the main function domain or a subdomain. Piecewise can describe the nature of the function rather than just be a way of expressing it. In addition, the term piecewise can also refer to any piecewise-defined function that applies to each piece but may not be applicable for the whole function. The piecewise or continuous differentiable behavior of a function is that each piece is differentiable within its domain, regardless of whether the whole function is differentiable at the intersection. When functions are analyzed convexly, the notation of a derivative can be replaced by that of the subderivative [17]. A function can indeed be described as “piecewise linear” or “piecewise continuous” or “piecewise differentiable” even if the pieces are not intervals. Equation (8) shows a piecewise notation that expresses the absolute value function as follows:(8)x=−x,ifx<0,+x,ifx≥0.To turn negative values into positive, the first function -x is used for allx values, and this invalidates the sign of the input value. The second function x is applied for all x values which could be larger than or equal to zero.


### 4.1. Linear Piecewise

The piecewise linear function is a collection of quadratic functions, each defined on a set of intervals of real numbers. A finite collection of such intervals must exist if the domain of the function will either need to be finite or be locally finite [20–23].There are several contexts in which piecewise linear functions are relevant. In general, the piecewise linear function can be defined on any vector space or affine space and simplicial complexes [18, 19, 24]. A linear function is not simply a linear transformation, but in these contexts, it is an affine linear function. There are two major subclasses of a piecewise linear function, called continuous and convex. As we can see from (9), every n-dimensional continuous function, in general, will have the following formulas:(9)Π∫PPℝn+1,flx=min∑∫Πmaxa,b∫Σa⟶⋅x⟶+b,Σ∫Pℝn+1,fx⟶=maxa⟶⋅x⟶+ba⟶,b∫Σ.


## 4.1. Linear Piecewise

The piecewise linear function is a collection of quadratic functions, each defined on a set of intervals of real numbers. A finite collection of such intervals must exist if the domain of the function will either need to be finite or be locally finite [20–23].There are several contexts in which piecewise linear functions are relevant. In general, the piecewise linear function can be defined on any vector space or affine space and simplicial complexes [18, 19, 24]. A linear function is not simply a linear transformation, but in these contexts, it is an affine linear function. There are two major subclasses of a piecewise linear function, called continuous and convex. As we can see from (9), every n-dimensional continuous function, in general, will have the following formulas:(9)Π∫PPℝn+1,flx=min∑∫Πmaxa,b∫Σa⟶⋅x⟶+b,Σ∫Pℝn+1,fx⟶=maxa⟶⋅x⟶+ba⟶,b∫Σ.


## 5. Conclusion

We provided a class of numerical approaches for approximating the solutions of differential equations in this study which are based on a modified version. When applied to ODEs, these approaches yield a continuous approximate solution that is accurate to order 2n at the nodes and to order n + 1 evenly over the interval. To extend the approaches to delay differential equations, some a posteriori tweaks boost the uniform accuracy to order 2n. Real-world issues are subjected to numerical testing and compared with other methodologies. Four methods have been presented for linear system integration of ODEs, and those methods require less memory and less computational effort than the low storage methods. The notations and formulations have been studied in such a way that subsequent development of the theory of differential equations will follow the structure of ODEs. We also presented that when a nonstandard piecewise function is used correctly, it produces fractional-order differential equations that converge to the original equation as the parameter leads to zero. Even though it was not developed from quadrature, the method has been demonstrated to be equivalent to quadrated-based methods. Nonuniform time steps can be used to accomplish the approach. The difference equation may appropriately describe the dynamics of the underlying fractional differential equation, as demonstrated in the formulas.



---

*Source: 4651084-2022-08-21.xml*