Mathematical Modeling of Memristors

Abstract

The memristor has quite a reputation as a missing circuit element. It is a powerful candidate for next-generation applications after being first implemented in HP’s laboratories. At this point, mathematical models were needed for the analysis of the memristor, and a lot of studies were done on this subject. In this chapter, mathematical modeling and simulations of the memristor device have been emphasized. Firstly, linear drift and nonlinear drift models have been described on the basic HP model. The window functions used in the nonlinear drift model have been widely examined. Different from HP model, the Simmons tunnel barrier and the threshold adaptive memristor model (TEAM) have been also mentioned. As a result, the most widely used modeling techniques have been described in detail.

1. Introduction

In the circuit theory, it refers to the existence of three basic circuit elements that define connections between basic circuit parameters such as current (i), voltage (v), charge (q), and magnetic flux (φ). These are resistor, inductor, and capacitor. However, a circuit element that determines the relationship between the charge and the magnetic flux is not defined. The fourth fundamental circuit element representing this relation was firstly presented by Chua in mathematical terms in 1971 with the name of the memristor (memory + resistor) [1]. In 2008, a group of researcher from HP laboratories announced that they were physically producing memristor [2]. In Figure 1, the relationship between fundamental circuit elements and basic circuit parameters is shown.

Figure 1.

Linkage between four fundamental circuit elements and basic circuit parameters such as current (i), voltage (v), charge (q), and magnetic flux (φ).

The relationship between current, voltage, charge, and flux for the memristor is given by:

$$\mathrm{v}\left(\mathrm{t}\right)=\mathrm{M}\left(\mathrm{q}\left(\mathrm{t}\right)\right)\mathrm{i}\left(\mathrm{t}\right)$$E1

$$\mathrm{M}\left(\mathrm{q}\right)=\mathrm{d\phi }\left(\mathrm{q}\right)/\mathrm{dq}$$E2

$$\mathrm{i}\left(\mathrm{t}\right)=\mathrm{W}\left(\mathrm{\phi }\left(\mathrm{t}\right)\right)\mathrm{v}\left(\mathrm{t}\right)$$E3

$$\mathrm{W}\left(\mathrm{\phi }\right)=\mathrm{dq}\left(\mathrm{\phi }\right)/\mathrm{d\phi }$$E4

where M(q) has the unity of resistance and W(φ) has the unity of conductance [1].

Memristor has properties that are different from other fundamental circuit elements and can only be seen in a memristor such as nonvolatile memory effect, passivity, and pinched hysteresis loop.

When Eqs. (1) and (2) are opened, Eqs. (5) and (6) can be written as follows:

$$\mathrm{v}\left(\mathrm{t}\right)=\mathrm{M}\left(\underset{-\infty }{\overset{\mathrm{t}}{\int }}\mathrm{i}\left(\mathrm{t}\right)\mathrm{dt}\right)\mathrm{i}\left(\mathrm{t}\right)$$E5

$$\mathrm{i}\left(\mathrm{t}\right)=\mathrm{W}\left(\underset{-\infty }{\overset{\mathrm{t}}{\int }}\mathrm{v}\left(\mathrm{t}\right)\mathrm{dt}\right)\mathrm{v}\left(\mathrm{t}\right)$$E6

Eqs. (5) and (6) show that the memristance value is related to the history of the current passing through the memristor. That is, when the current passing through the memristor is cut off, it remains at the value of the memristance value. The memristance value starts to change from the last value when it provides the current again to memristor. In other words, the memristor has a nonvolatile memory effect. On the other hand, memristor is not an element that stores energy [3, 4].

Memristor is similar to resistor with memory. It shows a nonlinear resistance characteristic that the charge parameter is state variable [5].

Another distinguishing feature of the memristor is that the I-V change shows the pinched hysteresis loop characteristic. A memristor fed by a bipolar periodic signal always exhibits a pinched hysteresis I-V characteristic that passes through the origin. As the frequency of the excitation signal increases, the hysteresis lobe area decreases monotonically. When the frequency tends to infinity, the pinched hysteresis loop shrinks toward a single-valued function [6, 7].

The memristor, with being a passive circuit element, has a unique ability to remember the state of resistance that it possesses by maintaining the relationship between voltage and current time integrals. Due to these features, they are being nominated for many different applications such as resistive memories, soft computing, neurocomputing, etc.

Different materials and techniques are used at the point of producing the memristor. The memristor structures produced from different materials can be given as example, such as titanium dioxide (TiO₂) memristor [2], zinc oxide memristor [8, 9], silicon oxide memristor [10], and GST (Ge₂Sb₂Te₅) memristor [11].

2. Modeling of memristor

2.1. HP memristor model

Memristor-based applications require a suitable model for analysis and simulation of the system. When looking at the literature, the HP memristor model where the memristor mechanism based on the drift of oxygen vacancies is widely used. The memristor model developed by HP Lab is composed of Pt/TiO₂/Pt structure as shown in Figure 2. Here, the TiO₂ layer in which the one side doped with positive charged-rich oxygen vacancies (TiO_2−x) is placed between two platinum layers [2].

Figure 2.

Structure of memristor reported by HP Lab [2].

The doped part of the TiO₂ layer exhibits a low resistance behavior, while the undoped part exhibits a high resistance behavior. As a result of the appropriate excitation on this structure, the ionic drift between the doped part and the undoped part results in a dynamic change in the width of the doped region. That is, the width of the doped region is taken as a state variable. As the width of the doped region approaches zero (w→0), the memristor goes to a high resistance state (HRS), and as the width of the doped region approaches D (w→D), the memristor goes to a low resistance state (LRS) as shown in Figure 3 [3].

Since the memristor’s dimensions are very small (in few nm), it causes a change in the doped region even with a small stimulation applied. Thus the resistance of the memristor varies between HRS and LRS [3].

2.2. Linear drift model

In the model known as the linear drift model, the relation between the current and the voltage of the memristor is defined by the following equation:

$$\mathrm{v}\left(\mathrm{t}\right)=\left[{\mathrm{R}}_{\mathrm{on}}\mathrm{x}\left(\mathrm{t}\right)+{\mathrm{R}}_{\mathrm{off}}\left(1-\mathrm{x}\left(\mathrm{t}\right)\right)\right]\mathrm{i}\left(\mathrm{t}\right)$$E7

$$\mathrm{x}\left(\mathrm{t}\right)=\frac{\mathrm{w}\left(\mathrm{t}\right)}{\mathrm{D}}\in \left(0,1\right)$$E8

where R_on and R_off are the values of the resistance for w(t) = D and w(t) = 0, respectively [2, 8, 9].

From Eq. (7), the value of memristance can be expressed by

$$\mathrm{M}\left(\mathrm{q}\left(\mathrm{t}\right)\right)=\frac{\mathrm{v}\left(\mathrm{t}\right)}{\mathrm{i}\left(\mathrm{t}\right)}={\mathrm{R}}_{\mathrm{on}}\mathrm{x}\left(\mathrm{t}\right)+{\mathrm{R}}_{\mathrm{off}}\left(1-\mathrm{x}\left(\mathrm{t}\right)\right)$$E9

As shown in Figure 3, the state of change of the memristor resistance is represented by x(t) value in Eq. (8). The speed of movement of the boundary between the doped layer and undoped layer is expressed as dx/dt with Eq. (10) [12, 13]:

$$\frac{\mathrm{dx}\left(\mathrm{t}\right)}{\mathrm{dt}}={\mathrm{\mu }}_{\mathrm{v}}\frac{{\mathrm{R}}_{\mathrm{on}}}{{\mathrm{D}}^2}\mathrm{i}\left(\mathrm{t}\right)$$E10

where μ_v is the average drift mobility of the charges. If Eq. (10) is taken integral for time, the following expression is derived:

$$\mathrm{x}\left(\mathrm{t}\right)={\mathrm{\mu }}_{\mathrm{v}}\frac{{\mathrm{R}}_{\mathrm{on}}}{{\mathrm{D}}^2}\mathrm{q}\left(\mathrm{t}\right)$$E11

Figure 3.

Representation of the HRS and LRS states of the memristor. (a) LRS and (b) HRS.

If Eq. (11) is put in its place in Eq. (9), the following memristance expression is achieved:

$$\mathrm{M}\left(\mathrm{q}\left(\mathrm{t}\right)\right)={\mathrm{R}}_{\mathrm{on}}{\mathrm{\mu }}_{\mathrm{v}}\frac{{\mathrm{R}}_{\mathrm{on}}}{{\mathrm{D}}^2}\mathrm{q}\left(\mathrm{t}\right)+{\mathrm{R}}_{\mathrm{off}}\left(1-{\mathrm{\mu }}_{\mathrm{v}}\frac{{\mathrm{R}}_{\mathrm{on}}}{{\mathrm{D}}^2}\mathrm{q}\left(\mathrm{t}\right)\right)$$E12

This expression can be written as

$$\mathrm{M}\left(\mathrm{q}\left(\mathrm{t}\right)\right)={\mathrm{R}}_{\mathrm{off}}\left(1-{\mathrm{\mu }}_{\mathrm{v}}\frac{{\mathrm{R}}_{\mathrm{on}}}{{\mathrm{D}}^2}\mathrm{q}\left(\mathrm{t}\right)\right)$$E13

if the left side of the total expression is neglected because R_on < < R_off [2].

Through this model, the characteristics of the memristor can be observed by the simulations as shown in the ongoing part. The following values are used for the simulations performed in this section. μ_v = 10⁻¹⁴ m²s⁻¹ V⁻¹, D = 10 nm, initial value of w is 3 nm, input signal V_input = V_o.sin(ωt) where V_o = 1 V and f = 1 Hz (ω = 2πf), R_on = 100 Ω, and R_off = 160 kΩ. Figure 4 shows the change of the current and voltage of the memristor with time for the given parameter values. The pinched hysteresis loop in the I-V plane shown in Figure 5 is one of the fingerprint characteristics of the memristor. Figure 6 shows the relationship between state variable and memristance. This indicates that the memristance depends on the state variable x. This figure also shows that the state variable is limited between 0 and 1. In Figure 7, the change of memristance with applied voltage is seen.

Figure 4.

Change of the current and voltage of the memristor with respect to time.

Figure 5.

I-V pinched hysteresis loop of the memristor.

Figure 6.

Change of the state variable and memristance of the memristor with respect to time.

Figure 7.

Change of the memristance of the memristor with respect to voltage.

As shown in Figure 8, as the frequency increases, the I-V pinched hysteresis loops become narrower. As the frequency increases toward infinity, the I-V characteristic seems to be a linear resistance characteristic. Figure 9 shows the variation of the I-V characteristic for different amplitude values of the excitation signal.

Figure 8.

I-V pinched hysteresis loops of the memristor for Vo = 1 V and different frequency values.

Figure 9.

I-V pinched hysteresis loops of the memristor for f = 1 Hz and different V_o values.

2.3. Nonlinear drift model and window functions

The linear drift model supposes that the state variable (x) of the memristor is proportional to the charge flowing through the memristor. This proportion is acceptable to the interface between the electrodes and the interface between the doped and undoped parts of the memristor. The position of the doped part changes with the applied input signal. Furthermore, the linear drift model assumes that the vacancies have the freedom to move along the all length of the memristor. These assumptions made in the HP model have been greatly simplified, neglecting some basic laws. The reported literature shows that the drift of vacancies is not linear in the region near the boundary interfaces. The reason is that even a small excitation signal can create a large electric field causing nonlinear drift of the vacancies near the boundary interfaces in the memristor. Another problematic situation related to the linear drift model is that the state variable (x) never reaches zero, indicating that oxygen vacancies are not present in the memristor. Similarly, the doped region cannot cover the entire length of the memristor, because there will be no undoped part and the memristor will not work in this way [3, 14, 15].

In order to provide nonlinearity for the boundary problems mentioned above, functions called window function are introduced. This function is implemented by rearranging the expression Eq. (10) as

$$\frac{\mathrm{dx}\left(\mathrm{t}\right)}{\mathrm{dt}}={\mathrm{\mu }}_{\mathrm{v}}\frac{{\mathrm{R}}_{\mathrm{on}}}{{\mathrm{D}}^2}\mathrm{i}\left(\mathrm{t}\right)\mathrm{f}\left(\mathrm{x}\left(\mathrm{t}\right)\right)$$E14

The function f(x) should have zero at the limits of the memristor (x = 0 and x = 1) and maximum value at the middle of the memristor (x = 0.5) [11]. An effective window function should satisfy the following conditions for modeling of nonlinearity [16]:

• The function should take account of the boundary situation at the top and bottom electrodes of the memristor.

• The function should provide nonlinear drift across the entire active area of the memristor.

• The function should ensure linkage between the linear and nonlinear drift models.

• The function should be scalable in the interval of f_max(x) can be obtained such that 0 ≤ f_max(x) ≤ 1.

• The function should include the control parameter to set the model.

Many different window functions are proposed as a result of the studies carried out in order to provide these criteria.

2.3.1. Joglekar’s window function

Joglekar’s window function can be given as

$$\mathrm{f}\left(\mathrm{x}\right)=1-{\left(2\mathrm{x}-1\right)}^{2\mathrm{p}}$$E15

where p is the control parameter which changes the flatness of the f(x) curve around its maximum value at x = 0.5 and is a positive integer [17].

In Figure 10, the change of window function proposed by Joglekar for different p values is shown. The characteristic of this function is similar to the rectangular window function by increasing p value, and the nonlinear drift effect is reduced. The disadvantage of Joglekar’s window function is the cling situation of the state variable at the boundaries, and it is difficult to change the window function due to the zero value at both boundaries. That is, the nonlinear drift problem is solved, but the boundary lock is not taken into account. When memristor arrives in R_on or R_off terminal condition, this state will be maintained forever due to zero value taken from the window function [13, 14].

Figure 10.

Joglekar window function for different p values.

2.3.2. Biolek’s window function

Biolek has introduced a window function that provides a solution for model errors (the cling situation of the state variable at the boundaries) of Joglekar’s window function. Biolek’s window function is expressed as follows:

$$\mathrm{f}\left(\mathrm{x}\right)=1-{\left(\mathrm{x}-\mathrm{stp}\left(-\mathrm{i}\right)\right)}^{2\mathrm{p}}$$E16

$$\mathrm{stp}\left(\mathrm{i}\right)=\left\{\begin{array}{l}1,\mathrm{i}\ge 0\\ 0,\mathrm{i}<0\end{array}\right.$$E17

where p is positive integer and i is the memristor current [18].

Figure 11 shows the variation of the window function proposed by Biolek for different p values. The proposed window function by Biolek depends not only on the state variable but also on the current flow through the memristor. Thus, the problem of boundary lock is resolved. However, this window function does not include the scalability factor, so the maximum value of the window function cannot be set to a lower or greater value [13].

Figure 11.

Biolek’s window function for different p values.

2.3.3. Prodromakis’ window function

Prodromakis’ window function is

$$\mathrm{f}\left(\mathrm{x}\right)=\mathrm{j}\left(1-{\left[{\left(\mathrm{x}-0.5\right)}^2+0.75\right]}^{\mathrm{p}}\right)$$E18

where p and j are a positive real number [16]. In Figure 12, the change of Prodromakis’ window function for j = 1 and different p values is shown. Figure 13 shows the variation of the Prodromakis’ window function for p = 10 and different j values.

Figure 12.

Prodromakis’ window function for j = 1 and different p values.

Figure 13.

Prodromakis’ window function for p = 10 and different j values.

Prodromakis proposes a solution for the scalability problem in the aforementioned by the presented window function. Prodromakis’ window function provides a connection to the linear dopant drift model for sufficiently large values of p. However, the model built by Prodromakis still contains the problem of boundary lock [13].

2.3.4. Zha’s window function

This function has been introduced by Zha as a new window model so that boundary lock, scalability, and nonlinear effects can be met at the same time. Zha’s window function is expressed as follows:

$$\mathrm{f}\left(\mathrm{x}\right)=\mathrm{j}\left(1-{\left[0.25{\left(\mathrm{x}-\mathrm{stp}\left(-\mathrm{i}\right)\right)}^2+0.75\right]}^{\mathrm{p}}\right)$$E19

where stp(i) is given in Eq. (17) and p and j are positive real numbers [13]. Zha’s window function for j = 1 and different p values in Figure 14 is shown. Figure 15 shows the Zha’s window function for p = 10 and different j values.

Figure 14.

Zha’s window function for j = 1 and different p values.

Figure 15.

Zha’s window function for p = 10 and different j values.

2.3.5. Comparison of window functions

In this section, the simulation results have been given over the nonlinear drift model of the window functions given in the previous sections. Figure 16 shows the change of window functions according to state variable x for p = 10 and j = 1 values. In Figure 17, I-V characteristics have been plotted for the window functions in Figure 16. In Figure 18, the change of the memristance of the memristor with the applied voltage for the window functions of Figure 16 has been presented. For simulations using these window functions, μ_v = 10⁻¹⁴ m²s⁻¹ V⁻¹, D = 10 nm, initial value of w is 3.145 nm, input signal V_input = V_o.sin(ωt) where V_o = 1.2 V and f = 1 Hz (ω = 2πf), R_on = 100 Ω, and R_off = 160 kΩ.

Figure 16.

Different window functions for p = 10 and j = 1.

Figure 17.

Change of I-V pinched hysteresis loops of the memristor for different window functions.

Figure 18.

Change of the memristance of the memristor with respect to voltage for different window functions.

2.4. Exponential model

Even in the models described so far, the nonlinearity of the large electric field in the memristor is still not taken into consideration. In [19], an exponential model that accounts the effect has been presented. In this model, the relation between the current of the memristor and the voltage is defined as follows:

$$\mathrm{i}=\mathrm{x}{\left(\mathrm{t}\right)}^{\mathrm{n}}\mathrm{\beta }\text{sinh}\left(\mathrm{\alpha }\mathrm{v}\left(\mathrm{t}\right)\right)+\mathrm{\chi }\left(\text{exp}\left(\mathrm{\gamma }\mathrm{v}\left(\mathrm{t}\right)\right)-1\right)$$E20

where β, α, χ, and γ are experimental fitting parameters. How the state variable can influence the current is determined by the n parameter [12, 19]. According to Eq. (20), when the model is ON state, asymmetrical switching behavior is shown (sinh part). When the OFF state, the exponential part of the Eq. (20) has the dominant part of the current, which is similar to an ideal PN junction [12, 14].

In this model, the differential equation of state variable is written as

$$\frac{\mathrm{dx}\left(\mathrm{t}\right)}{\mathrm{dt}}=\mathrm{a}\cdot \mathrm{v}{\left(\mathrm{t}\right)}^{\mathrm{m}}\mathrm{f}\left(\mathrm{x}\right)$$E21

where a and m are fitting parameters. f(x) can be any window function [14].

2.5. Simmons tunnel barrier model

The models described so far were based on the HP model, which consisted of two regions, each of which was modeled as resistance. But Pickett presented another physical model of the memristor as an alternative to the HP model, consisting of a resistor and an electron tunnel barrier in series [20].

Figure 19 shows the memristor structure of the Simmons tunnel barrier model where w is the tunneling barrier and R_s is the channel resistance.

Figure 19.

Memristor structure of Simmons tunnel barrier model.

w is the state variable of the model and can be written as

$$\frac{\mathrm{dw}\left(\mathrm{t}\right)}{\mathrm{dt}}=\left\{\begin{array}{l}{\mathrm{f}}_{\mathrm{off}}\text{sinh}\left(\frac{\mathrm{i}}{{\mathrm{i}}_{\mathrm{off}}}\right)\text{exp}\left[-\text{exp}\left(\frac{\mathrm{w}-{\mathrm{a}}_{\mathrm{off}}}{{\mathrm{w}}_{\mathrm{c}}}-\frac{\left|\mathrm{i}\right|}{\mathrm{b}}\right)-\frac{\mathrm{w}}{{\mathrm{w}}_{\mathrm{c}}}\right],\mathrm{i}>0\\ {\mathrm{f}}_{\mathrm{on}}\text{sinh}\left(\frac{\mathrm{i}}{{\mathrm{i}}_{\mathrm{on}}}\right)\text{exp}\left[-\text{exp}\left(\frac{\mathrm{w}-{\mathrm{a}}_{\mathrm{on}}}{{\mathrm{w}}_{\mathrm{c}}}-\frac{\left|\mathrm{i}\right|}{\mathrm{b}}\right)-\frac{\mathrm{w}}{{\mathrm{w}}_{\mathrm{c}}}\right],\mathrm{i}<0\end{array}\right.$$E22

where f_off, f_on, a_off, a_on, i_off, i_on, and b are fitting parameters [20]. The f_on value has an amplitude order greater than f_off. It is also effective in changing w in both parameters. The i_on and i_off parameters effectively limit the current threshold. In this model, a window function is not required, because a_off and a_on values force upper and lower bounds of x, respectively. Although this model is the most accurate model for a memristor, it is a nongeneric model that is defined for a particular type of memristor, which has a nonobvious relationship between current and voltage [14].

2.6. ThrEshold Adaptive Memristor (TEAM) model

TEAM model is a memristor model with several assumptions for analysis simplification and computational efficiency. These assumptions are as follows:

1. There is no change in the status variable for values below a certain threshold value.

2. A polynomial relationship is established between the current of memristor and the internal state drift derivative instead of the exponential dependence.

Taking these assumptions into account, the derivation of the state variable is written as

$$\frac{\mathrm{dw}\left(\mathrm{t}\right)}{\mathrm{dt}}=\left\{\begin{array}{l}{\mathrm{k}}_{\mathrm{off}}\cdot {\left(\frac{\mathrm{i}\left(\mathrm{t}\right)}{{\mathrm{i}}_{\mathrm{off}}}-1\right)}^{{\mathrm{\alpha }}_{\mathrm{off}}}\cdot {\mathrm{f}}_{\mathrm{off}}\left(\mathrm{w}\right),0<{\mathrm{i}}_{\mathrm{off}}<\mathrm{i}\\ 0,{\mathrm{i}}_{\mathrm{on}}<\mathrm{i}<{\mathrm{i}}_{\mathrm{off}}\\ {\mathrm{k}}_{\mathrm{on}}\cdot {\left(\frac{\mathrm{i}\left(\mathrm{t}\right)}{{\mathrm{i}}_{\mathrm{on}}}-1\right)}^{{\mathrm{\alpha }}_{\mathrm{on}}}\cdot {\mathrm{f}}_{\mathrm{on}}\left(\mathrm{w}\right),\mathrm{i}<{\mathrm{i}}_{\mathrm{on}}<0\end{array}\right.$$E23

where k_off (k_off > 0,) k_on (k_on < 0),α_off, and α_on are constants, i_off and i_on are current thresholds, and w is the effective electric tunnel width. Dependency on state variable w by f_off(w) and f_on(w) functions is provided. These functions can be thought of as window functions to limit the state variable between w_on and w_off. If we assume that the memristance changes linearly with w as in Eq. (7), the relationship between current and voltage can be written as

$$\mathrm{v}\left(\mathrm{t}\right)=\left[{\mathrm{R}}_{\mathrm{on}}+\frac{{\mathrm{R}}_{\mathrm{off}}-{\mathrm{R}}_{\mathrm{on}}}{{\mathrm{w}}_{\mathrm{off}}-{\mathrm{w}}_{\mathrm{on}}}\left(\mathrm{w}-{\mathrm{w}}_{\mathrm{on}}\right)\right]\mathrm{i}\left(\mathrm{t}\right)$$E24

If we assume that the memristance changes exponentially with w, the relationship between current and voltage can be written as

$$\mathrm{v}\left(\mathrm{t}\right)={\mathrm{R}}_{\mathrm{on}}{\mathrm{e}}^{\left(\frac{\mathrm{\lambda }}{{\mathrm{w}}_{\mathrm{off}}-{\mathrm{w}}_{\mathrm{on}}}\left(\mathrm{w}-{\mathrm{w}}_{\mathrm{on}}\right)\right)}\cdot \mathrm{i}\left(\mathrm{t}\right)$$E25

where λ is fitting parameter [21]:

$$\mathrm{\lambda }=\text{ln}\left(\frac{{\mathrm{R}}_{\mathrm{off}}}{{\mathrm{R}}_{\mathrm{on}}}\right)$$E26

3. Conclusion

This chapter describes the mathematical modeling and simulation of the memristor device. Firstly, brief information about the historical development of the memristor has been given. The emergence of the memristor idea and the formation of mathematical theory have been mentioned. Then, information about the realization of the memristor as a physical element and the HP memristor model has been given.

In the memristor applications, the memristor device must be mathematically modeled correctly for analysis and simulation studies. For this reason, mathematical modeling and modeling methods of the memristor have been emphasized.

We mainly focus on five different models such as linear drift model, nonlinear drift model, exponential model, Simmons tunnel barrier model, and TEAM model. In addition, the different window functions proposed for the nonlinear drift model have been examined. We provided simulation results for some of the models reviewed. The effects on the I-V characteristics of the window functions have been shown graphically with simulation results.