Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3226 https://internationalpubls.com Numerical Accuracy of Euler-Maruyama and Milstein Methods for Solving Nonlinear Stochastic Differential Equation Eram Naaz1 and Najmuddin Ahmad2* 1,2Department of Mathematics & Statistics, Integral University, Lucknow- 226026, India. *Coressponding author e-mail: najmuddinahmad33@gmail.com 2 Author e-mail: eramphdiul@gmail.com 1 Article History: Received: 04/01/2025 Revised: 06/02/2025 Accepted: 11/03/2025 Abstract: In this paper we examine the numerical accuracy of two widely used methods-the Euler–Maruyama and Milstein methods for solving Stochastic Differential Equations (SDEs). In this study we deal with two nonlinear SDEs, whose exact solutions are derived using Itô’s calculus. Using uniform discretization conditions, the comparative study is made on the exact solutions, approximate solutions, and absolute errors, supported with graphs and tables. This study focuses on selecting reliable numerical methods for accurately solving nonlinear stochastic systems. Keywords: Stochastic Differential Equations, Itô Lemma, Euler- Maruyama Method, Milstein Method. AMS Subject Classification: 65A05, 65B15 1. Introduction While modeling complex real-world systems that have randomness and uncertainty, like stock markets, neurons, and ecosystems that change suddenly, the traditional deterministic models like ordinary differential equations are insufficient because of their inability to take into account stochastic fluctuations and unpredictable behaviors. This gap has been addressed by stochastic differential equations (SDEs) by capturing both the deterministic trends and the random noise that affects the system’s behavior [1]. For example, in quantitative finance, SDEs model the uncertain behavior of stock prices and interest rates better than the traditional models [1], [13]. Similarly, SDEs help to simulate the effects of random outside stimuli and mutations in neuroscience and epidemiology [1]. SDEs incorporate stochastic elements like Brownian motion, making them more realistic and detailed for the systems having randomness or uncertainty, making them more valuable across various disciplines like finance, physics, biology, engineering, and signal processing [1]. Fundamental research on stochastic processes was driven by Einstein's 1905 analysis of Brownian motion, which provided a probabilistic interpretation of diffusion phenomena and served as the theoretical foundation for stochastic differential equations (SDEs) [1]. Norbert Wiener made this perspective very formal by developing the Wiener process, which provided a mathematical model for stochastic perturbations with stationary and independent mailto:najmuddinahmad33@gmail.com mailto:eramphdiul@gmail.com Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3227 https://internationalpubls.com increments in continuous time [1]. Kiyosi Itô's contributions to stochastic calculus, including Itô's lemma and Itô integral, enabled the differential analysis of non-deterministic systems driven by Brownian motion and resulted in the modern formulation of SDEs [1]. Paul Lévy added to this theory by introducing Lévy processes, which can handle jump discontinuities and extend Brownian motion to a wider range of stochastic drivers [1]. These contributions together formed the theoretical basis of SDEs, which let us model systems that are affected by both jump-type and continuous randomness. In various scientific and engineering fields, SDEs are now essential for the analysis and simulation of complex systems [1], [13]. 2. Stochastic Differential Equation A crucial tool for simulating systems with inherent unpredictability is the stochastic differential equation (SDE). But they still face issues like analytical intractability, especially with nonlinear SDEs, where closed-form solutions are very difficult to find. Because of this limitation, traditional stochastic techniques are unable to fully capture the complex relationship between stochastic perturbations and deterministic dynamics that happens in the real-world [2]. The importance of advanced numerical techniques in solving complicated mathematical models is further illustrated by the promising results of semi-analytical approaches like the Laplace–Adomian Decomposition Method in handling nonlinear integral equations [10]. As a result of this, the numerical approximation of SDEs is now a crucial part of stochastic modeling. The Euler–Maruyama and Milstein methods are renowned among the developed numerical methods as fundamental first-order approaches [3], [7], [13]. The Euler- Maruyama method is widely used due to its conceptual simplicity and ease of use, as it is a natural stochastic analogue of the classical Euler method [13]. However, its accuracy in finer simulations is limited by its relatively low order of strong convergence [3]. This issue is solved by the Milstein method by adding an additional term involving the diffusion coefficient's derivative, improving the quality of the approximation, particularly in systems with multiplicative noise or non-trivial diffusion structures [3], [6], [16]. Modified numerical methods have been researched because classical methods like Euler-Maruyama and Milstein may not be able to hold important characteristics like positivity or stability in situations involving super-linear growth or delay effects [9], [14]. The conceptual basis for stochastic numerical analysis was established by Kloeden and Platen [13], Higham [3], Burrage [7], and others in their seminal works [3], [7], [13], which thoroughly investigated the theoretical foundations and convergence analyses of these methods. After these developments, even now most of the existing comparative literature focuses mainly on linear or simplified SDEs, frequently without access to known analytical solutions. This makes it very difficult to accurately measure numerical error and see how well the methods work in more complicated situations [4], [15]. Identifying this gap, the current study conducts a systematic comparison of the Euler–Maruyama and Milstein methods using two carefully selected nonlinear SDEs, each with precise closed-form solutions obtained via Itô calculus [1],[6]. These exact solutions are not only difficult to obtain but also provide strict criteria to determine the precision, stability, and convergence behavior of both approaches under uniform time steps. This work contributes to a better understanding of the performance and limitations of these widely used numerical methods by concentrating on systems where both stochastic effects Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3228 https://internationalpubls.com and nonlinearity are important. Thus, provide a better understanding of their applicability in real-world modeling situations [2],[13]. 2.1 Methods In this paper we consider the general form of a one-dimensional SDE with dX(t) = α(X(t), t) dt + β(X(t), t) dW(t), X(t0) = X0, t0 ≤ t ≤ T (1) where: X(t) is the unknown stochastic process. α(X(t), t) is the drift coefficient. β(X(t), t) is the diffusion coefficient. W(t) is the Wiener process (also called standard Brownian motion), representing the source of randomness. Robert Brown noticed the jittery motion of microscopic particles suspended in a fluid in 1827, which led to the development of the theory of Brownian motion. Later, Norbert Wiener's work gave this physical phenomenon a rigorous mathematical foundation by constructing the Wiener process, which served as the basis for modern stochastic calculus [1]. The Wiener process W(t) on [0,T] is characterized by three essential properties: i. W(0)=0, which means that it begins deterministically at zero (with probability 1). ii. For any 0 ≤ s < t ≤ T, the increment W(t)−W(s) is N(0, t-s) distributed. iii. The increments W(t) − W(s) and W(v) − W(u) are statistically independent for disjoint intervals [s,t] and [u,v] respectively, where 0 ≤ s < t < u < v ≤ T. Researchers like Xuerong Mao and Chenggui Yuan have made significant contributions to the study of stochastic systems. For instance, they have created theoretical frameworks for stochastic differential equations (SDEs) with Markovian switching and nonlinear dynamics [17] and also developed numerical methods that preserve structural properties in financial SDEs with superlinear coefficients [14]. Integrating equation (1) from t0 to t, we have the following integral form: X(t) = X(t0) + ∫t0 t α(s, X(s))ds + ∫t0 t β(s, X(s))dW(s) The first integral ∫t0 t α(s,X(s))ds is a standard (Riemann or Lebesgue) integral that represent the deterministic part (drift). The second integral ∫t0 t β(s, X(s))dW(s) is a stochastic (Itô) integral, that represents the random fluctuations introduced by the Brownian motion W(t). The last term in the integral formulation of a stochastic differential equation is known as the Itô integral; it is a key idea in stochastic calculus. The Itô formula is a stochastic variant of Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3229 https://internationalpubls.com the classical calculus chain rule; it is used to analytically examine such equations, especially when transforming functions of stochastic processes [1], [13]. We create a new process. U = f(t, Xt), where, f is a smooth function and the differential of U is given by: dU = (∂f/∂t)dt + (∂f/∂x)dXt + 1/2(∂2f/∂x2)(dXt)^2 Substituting the expression for dXt into this expansion and applying Itô’s product rules— Itô multiplication table We have, dU =[ (∂f/∂t) + α(∂f/∂x) + (1/2)β2(∂2f/∂x2)]dt + β(∂f/∂x)dWt This formulation is the basis for Itô–Taylor expansions, which are fundamental in developing high-precision numerical approaches for solving stochastic differential equations [3], [13]. Over the years, various studies have been done comparing the effectiveness of numerical techniques in solving ordinary differential equations with initial value problems [11]. These basic ideas continue to influence the development and assessment of more advanced methods such as Euler–Maruyama and Milstein for stochastic systems. 2.2 Euler Maruyama Method A simple technique to get approximate answers to stochastic differential equations is the Euler–Maruyama method. It is like a traditional Euler method for solving deterministic ordinary differential equations (ODEs) [3], [13]. This method extends the concept of Taylor expansion to incorporate stochastic processes, forming the basis for simulating Itô-type SDEs in real-world situations where analytical solutions are hard or impossible to find [13]. As explained in [13] the scheme is obtained by truncating the Itô-Taylor expansion after the first- order terms (drift and diffusion). Let us consider the Itô stochastic differential equation of the form: dXt = α(Xt, t) dt + β(Xt, t)dWt, Xt0 = X0, . dt dWt dt 0 0 dWt 0 dt Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3230 https://internationalpubls.com where α(Xt, t) is the drift coefficient, β(Xt, t) is the diffusion coefficient, and Wt is the Wiener process or Brownian motion. To numerically approximate this SDE, we divide the time interval [t0, T] into N equal steps of size Δt. We define the Euler–Maruyama iteration scheme using the formula: Yn+1 = Yn + α(Yn, tn)Δt + β(Yn,tn)ΔWn , where ΔWn = Wtn+1 − Wtn is the Brownian increment over the interval, and E[ΔWn] = 0 and Var(ΔWn) = Δt. This scheme is derived by reducing the stochastic Taylor expansion at the lowest-order terms [2]. The deterministic part, α(tn,Yn)Δt, is derived from the standard Euler method, while the stochastic term β(tn,Yn)ΔWn accounts for the random fluctuations. The Euler–Maruyama method achieves a strong convergence of order 0.5, as evidenced by the decreasing expected error between the exact and numerical solutions with √Δt. On the other hand, as discussed in [3], [13], its weak convergence (for expectation errors) is of order 1. When the diffusion term β ≡ 0, the method reduces to the deterministic Euler scheme with a stronger convergence rate of order 1 [1], [3], [13]. Despite its limitations, the Euler–Maruyama method remains popular for modeling systems with randomness. Because of its simplicity and ease of use, it is widely used in finance, physics, biology, and other applied areas [1], [8], [13]. 2.3 Milstein Method The Milstein method is an advanced numerical method developed to improve upon the Euler–Maruyama method by adding an additional correction term from the second-order Itô- Taylor expansion for approximating solutions of stochastic differential equations (SDEs), hence achieving a strong convergence order of 1 [3], [6], [13]. The Euler–Maruyama method approximates the SDE dXt = α(Xt,t) dt + β(Xt,t) dWt by discretizing the drift and diffusion terms, while the Milstein method improves this by adding the derivative of the diffusion term. The Milstein approximation for this SDE, over a uniform discretization tn = nΔt, is given by: Xn+1 = Xn + α(Xn)Δt + β(Xn)ΔWn + (1/2)β(Xn)β'(Xn)[(ΔWn)2 − Δt], where ΔWn = Wtn+1 − Wtn and β′(Xn) is the partial derivative of β with respect to the spatial variable X. The additional term, (1/2)β(Xn)β'(Xn)[(ΔWn)2 − Δt], captures the random variations of Brownian motion more accurately by including the Itô–Stratonovich adjustment that the Euler–Maruyama method ignores. When dealing with nonlinear stochastic systems (where β depends on Xt), the Milstein method performs significantly better than Euler–Maruyama [6], [13], [16]. When the noise is additive (β' = 0), both the methods give the same result. Due to its ability to capture higher-order stochastic dynamics, it is considered a better option in applications where fine-grained statistical simulation, strong convergence, and path accuracy Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3231 https://internationalpubls.com are essential. This makes it valuable in fields like quantitative finance, stochastic physics, and computational biology, where capturing complex stochastic behavior is essential [6], [13]. 3. Numerical Implementation and Comparative Analysis For the numerical implementation and comparative analysis, we consider two distinct nonlinear stochastic differential equations (SDEs) with known exact solutions. This is done to evaluate the accuracy and convergence behavior of the Euler–Maruyama and Milstein methods when applied to systems with nonlinear drift and diffusion structures [4], [15]. The selected SDEs are: Example 1: dX(t) = aXt (1 - Xt 2)dt + b(1 - Xt 2)dWt where, a = 1, b = 1.5, T = 1, X0 = 0.5, t ϵ [0, T] The exact solution given by: X(t) = tanh(tanh⁻¹(X0) + (a + b²)t + bWt) Example 2: dX(t) = [(2/5)Xt^(3/5) + 5Xt^(4/5)]dt + Xt^(4/5)dWt where, T = 1, X0 = 1, t ϵ [0, T] The exact solution given by: X(t) = (1 + t + (1/5)Wt) [8] Euler–Maruyama and Milstein schemes are used to numerically integrate both SDEs over a uniform time grid. The computed results are compared with the exact solution at various time points, and the results are compiled in tabular form to quantitatively evaluate the performance of both methods. This comparison enables a deeper understanding of the methods’ accuracy in handling nonlinear dynamics and their responsiveness to the structure of the diffusion term. Along with the tabulated data, graphical comparisons are presented to show the exact solution and the trajectories generated by each numerical method. These plots show how accurately the Milstein and Euler–Maruyama approximations follow the exact solution over time. By including both individual and combined trajectories, the figures effectively illustrate variations in pathwise behavior and provide insightful information about the relative accuracy of each scheme. Previous research on ordinary differential equations has shown that increasing the algorithmic order of fifth-order Runge-Kutta methods can enhance the accuracy of numerical results [12]. These classical results provide conceptual understanding of numerical analysis. Due to the non-differentiable nature of Brownian motion, the stochastic differential equations (SDEs) present fundamentally different challenges. As a result, even the most advanced methods, like Euler Maruyama and Milstein achieve significantly lower convergence orders Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3232 https://internationalpubls.com than their ODE counterparts [3], [13]. The result of ODEs [12] remain valuable for understanding and comparing how different numerical methods perform in stochastic situations. The table below shows the computed solutions for Example 1 and Example 2, obtained using MATLAB. Table 1 - Comparative Evaluation of Euler–Maruyama and Milstein Methods for Example 1 Figure 1 - Exact vs Euler Maruyama vs Milstein Approximation Time Exact Solution Euler Maruyama Milstein 0.0020 0.523426969817 0.520918917534 0.521053465775 0.1020 0.548354657801 0.390362133547 0.392628655112 0.2020 0.667506138063 0.356521561587 0.361200653554 0.3020 0.945028337189 0.832924520116 0.847678837873 0.4020 0.964242605866 0.876321176132 0.891026800035 0.5020 0.983133376703 0.948884254395 0.944217419428 0.6020 0.978290012959 0.928645239899 0.926191258282 0.7020 0.988486234074 0.960326979129 0.958629463677 0.8020 0.998297123373 0.994057619480 0.993773765945 0.9020 0.998777951262 0.995658925946 0.995504785803 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3233 https://internationalpubls.com Table 2 - A Comparative Error Study of Euler–Maruyama and Milstein Methods for Example 1 The mean absolute error for Euler Maruyama is 0.079092207934 and mean absolute error for Milstein is 0.076363830162. Figure 2 - Absolute Error Comparison Time Exact Solution Absolute Error (Euler M.) Absolute Error (Milstein) 0.0020 0.523426969817 0.002508052283 0.002373504042 0.1020 0.548354657801 0.157992524254 0.155726002689 0.2020 0.667506138063 0.310984576476 0.306305484509 0.3020 0.945028337189 0.112103817073 0.097349499316 0.4020 0.964242605866 0.087921429734 0.073215805831 0.5020 0.983133376703 0.034249122308 0.038915957275 0.6020 0.978290012959 0.049644773060 0.052098754677 0.7020 0.988486234074 0.028159254945 0.029856770397 0.8020 0.998297123373 0.004239503893 0.004523357428 0.9020 0.998777951262 0.003119025316 0.003273165459 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3234 https://internationalpubls.com Table 3 - Comparative Evaluation of Euler–Maruyama and Milstein Methods for Example 2 Figure 3 - Exact vs Euler Maruyama vs Milstein Approximation Time Exact Solution Euler Maruyama Milstein 0.0020 0.974885356941 0.974862988485 0.974897937607 0.1020 1.471179176072 1.470652743818 1.471123994331 0.2020 3.990556304959 3.977430661365 3.989486970767 0.3020 4.799183231921 4.783163300458 4.797980760965 0.4020 7.377947326403 7.357715839617 7.375795333862 0.5020 9.418679276993 9.397619125580 9.415924333430 0.6020 15.557996089632 15.535273112617 15.552684640900 0.7020 13.898007189249 13.863297287643 13.894091286272 0.8020 17.085661386138 17.051310284758 17.080860625168 0.9020 14.892091213402 14.859788368000 14.888848395348 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3235 https://internationalpubls.com Table 4 - A Comparative Error Study of Euler–Maruyama and Milstein Methods for SDE 2 The mean absolute error for Euler Maruyama is 0.019507283937 and mean absolute error for Milstein is 0.002451743439 . Figure 4 - Absolute Error Comparison Time Exact Solution Absolute Error (Euler M.) Absolute Error (Milstein) 0.0020 0.974885356941 0.000022368456 0.000012580666 0.1020 1.471179176072 0.000526432254 0.000055181741 0.2020 3.990556304959 0.013125643594 0.001069334192 0.3020 4.799183231921 0.016019931463 0.001202470956 0.4020 7.377947326403 0.020231486786 0.002151992541 0.5020 9.418679276993 0.021060151413 0.002754943563 0.6020 15.557996089632 0.022722977015 0.005311448732 0.7020 13.898007189249 0.034709901606 0.003915902977 0.8020 17.085661386138 0.034351101380 0.004800760970 0.9020 14.892091213402 0.032302845402 0.003242818054 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3236 https://internationalpubls.com 4._Conclusion In this study, we compare the Euler–Maruyama and Milstein methods using two nonlinear Stochastic Differential Equations (SDEs) with known exact solutions. Overall, the results show that the Milstein method is usually more accurate, but its performance depends largely on the SDE's structural properties, especially the diffusion term's nonlinearity and differentiability. Milstein clearly gives better results than Euler–Maruyama when the diffusion term is regular and highly nonlinear. The numerical experiments also demonstrate that Milstein is more accurate in both cases, but the amount of improvement is different. Thus, the structural properties of the SDE specifically the forms of the drift and diffusion coefficients, should be considered when selecting a suitable numerical method. Acknowledgement: IU/R&D/2025-MCN0003831 office of research and development Integral University, Lucknow. References [1] B. Øksendal, Stochastic Differential Equations: An Introduction with Applications, 6th ed. Berlin, Germany: Springer, 2003. [2] B. J. Akinbo, T. Faniran, and E. O. Ayoola, “Numerical solution of stochastic differential equations,” Int. J. Adv. Res. Sci. Eng. Technol., vol. 2, no. 5, pp. 608–616, 2015. [3] D. J. Higham, “An algorithmic introduction to numerical simulation of stochastic differential equations,” SIAM Rev., vol. 43, no. 3, pp. 525–546, 2001. [4] E. Buckwar and T. Sickenberger, “A comparative linear mean-square stability analysis of Maruyama- and Milstein-type methods,” Math. Comput. Simul., vol. 81, no. 6, pp. 1110–1127, 2011. [5] E. S. Rahman, “Euler-Maruyama’s method for numerical solution of Ornstein– Uhlenbeck’s equation,” Dept. Math., Fac. Math. Nat. Sci., Padang State Univ., 2020. [6] G. N. Mil’shtein, “Approximate integration of stochastic differential equations,” Theory Probab. Appl., vol. 19, no. 3, pp. 557–562, 1974. [7] K. Burrage, P. Burrage, and T. Mitsui, “Numerical solutions of stochastic differential equations—Implementation and stability issues,” J. Comput. Appl. Math., vol. 125, pp. 171– 182, 2000. [8] M. Bayram, T. Partal, and G. O. Buyukoz, “Numerical methods for simulation of stochastic differential equations,” Math. Comput. Simul., 2018. [9] M. T. Griggs, K. Burrage, and P. M. Burrage, “Magnus methods for stochastic delay- differential equations,” arXiv preprint arXiv:2506.16908, Jun. 2025. [10] B. M. Singh and N. Ahmad, “Numerical solution of Volterra nonlinear integral equation by using Laplace Adomian decomposition method,” Int. J. Appl. Math., vol. 35, no. 1, pp. 39–48, Mar. 2022. [11] S. Charan and N. Ahmad, “A comparative study on numerical solution of ordinary differential equation by different method with initial value problem,” Int. J. Recent Sci. Res., vol. 8, no. 10, pp. 21134–21139, Oct. 2017. [12] N. Ahamad and S. Charan, “Study of numerical solution of fourth order ordinary differential equations by fifth order Runge-Kutta method,” Int. J. Sci. Res. Sci. Eng. Technol. (IJSRSET), vol. 6, no. 1, pp. 230–238, Jan.–Feb. 2019. [13] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations. Berlin, Germany: Springer, 1992. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 32 No. 10s (2025) 3237 https://internationalpubls.com [14] S. Deng, C. Fei, W. Fei, and X. Mao, “Positivity-preserving truncated Euler and Milstein methods for financial SDEs with super-linear coefficients,” arXiv:2410.05614 [math.NA], Oct. 8, 2024. [15] S. E. Fadugba, B. J. Adegboyegun, and O. T. Ogunbiyi, “On convergence of Euler- Maruyama and Milstein scheme for solution of stochastic differential equations,” Int. J. Appl. Math. Model., vol. 1, no. 1, pp. 9–15, 2013. [16] S. J. Kayode and A. A. Ganiyu, “Effect of varying step sizes in numerical approximation of stochastic differential equations using one step Milstein method,” Appl. Comput. Math., vol. 4, no. 5, pp. 351–362, 2015. [17] X. Mao and C. Yuan, Stochastic Differential Equations with Markovian Switching. London, U.K.: Imperial College Press, 2006.