This study evaluates the accuracy and stability of the Euler and Runge-Kutta (RK4) methods using Wolfram Mathematica. The algorithms were tested on three models: a linear baseline, a stiff nonlinear system, and a real-world atmospheric free-fall model based on NASA Glenn Research Center physics. Results show that while Euler’s method is suitable for simple linear cases, it incurs significant drift in nonlinear systems, such as Exercise 2 where it produced a 39.2% error. In contrast, the RK4 method maintained high accuracy all through the three exercises. The findings confirm that RK4’s ) convergence is essential for reliable engineering and aerodynamic simulations where nonlinear drag forces are present.
O(
4
h
Keywords: Euler’s Method, Runge-Kutta Method (RK4), Ordinary Differential Equations (ODEs), Rule-Based Programming.
Nomenclature
Nomenclature
y’ or v’:The rate of change (Derivative).
h:The step size (Temporal resolution).
:Thenumericalapproximationatstepn.
h:The step size (Temporal resolution).
y
n
Introduction
Introduction
A differential equation is an equation that involves an unknown function together with one or more of its derivatives with respect to one or more independent variables. When the unknown function depends on a single independent variable, the equation is called an ordinary differential equation (ODE). When the unknown function depends on two or more independent variables and the derivatives involved are partial derivatives , the equation is called a partial differential equation (PDE).
∂
There are three main goals of studying differential equations ref[1,2]:
◼
To formulate mathematical models of real-life phenomena involving change such as physical, biological, and engineering systems.
◼
To analyze and solve these equations in order to understand the behavior of the modeled systems.
◼
To predict and control system behavior by interpreting solutions and using them to make decisions.
The following are examples of differential equations, ordinary and partial differential equations respectively.
dy
dx
∂u(x,t)
∂t
2
∂
∂
2
x
For relatively simple equations, a variety of analytical methods exist which includes; method of separation of variables, integrating factors, Bernoulli’s equation, and exact differential equations. These techniques exploit the structural form of an equation to obtain explicit closed-form solutions. However, many differential equations arising in practical applications are nonlinear and mathematically complex, making exact analytical solutions difficult or impossible to obtain.
As a result these reasons, numerical methods and semi-analytical techniques become necessary for obtaining meaningful and practical solutions. Among the most fundamental numerical methods are Euler’s method and Runge-Kutta methods which approximate solutions by advancing the system in discrete steps.
The primary goal of this work is to develop a rule-based implementation of Euler’s method and Runge–Kutta methods for solving ordinary differential equations, along with computation of exact solutions to obtain absolute errors, using the Wolfram Language, with emphasis on their algorithmic structure and implementation.
As a result these reasons, numerical methods and semi-analytical techniques become necessary for obtaining meaningful and practical solutions. Among the most fundamental numerical methods are Euler’s method and Runge-Kutta methods which approximate solutions by advancing the system in discrete steps.
The primary goal of this work is to develop a rule-based implementation of Euler’s method and Runge–Kutta methods for solving ordinary differential equations, along with computation of exact solutions to obtain absolute errors, using the Wolfram Language, with emphasis on their algorithmic structure and implementation.
Objectives of the work
Objectives of the work
The specific objectives of the work are to:
◼
Formulate Euler and Runge–Kutta methods to solve selected ODEs in an explicit step-by-step computational framework.
◼
implement the methods using concise Wolfram Language code.
◼
compare the accuracy and stability of Euler’s Method (First-Order) and the Classical Runge-Kutta Method (RK4) when solving linear, nonlinear, and real-world physical systems.
Novelty of the work
Novelty of the work
The novelty of this work lies in presenting classical numerical methods using rule-based computational steps. This approach highlights how simple iterative rules can be applied to approximate solutions of differential equations, making the methods transparent, reproducible, and suitable for experimentation in the Wolfram Language.
Euler method
Euler method
Let us consider a first-order ODE.
dy
dx
(
1
)subject to initial condition .
The numerical methods for solving differential equations are generally obtained initially by partitioning the independent variable into finite subintervals at steps ,,,...,, having equal step size for .
In terms of Taylor’s series expansion, the solution function of Equation (1) is obtained using
y()=
x
0
y
0
The numerical methods for solving differential equations are generally obtained initially by partitioning the independent variable
x∈[,]
x
0
x
n
n
x
0
x
1
x
2
x
n-1
x
n
h=Δx=-
x
i
x
i-1
i=1,2,...,n-1
In terms of Taylor’s series expansion, the solution function
y(x)
y(x+Δx)=y(x)+Δx+y+y+...
dy
dx
2
Δx
2!
2
d
d
2
x
3
Δx
3!
3
d
d
3
x
(
2
)y(x+h)=y(x)+h+y+y+...
dy
dx
2
h
2!
2
d
d
2
x
3
h
3!
3
d
d
3
x
(
3
)Equation (3) reduces to Euler’s equation if we retain the first-order approximation and neglect higher-order terms giving
y(x+h)=y(x)+h
dy
dx
(
4
)substituting =F(x,y) from Equation (1) gives;
dy
dx
y(x+h)=y(x)+hF(x,y)
(
5
) For discrete steps , , ..., Equation (5) may be written as
x
0
x
1
x
n-1
y
i+1
y
i
x
i
y
i
(
6
) For . Here Equation (6) is known as the Euler method or Euler-Cauchy method.
i=0,1,...,n-1
Algorithm: Euler method
Algorithm: Euler method
Input: Differential equation, with initial condition ; Domain, ; Step size, ; Number of subintervals, . Step(i): Discretization of the domain Define a uniform grid of points =+i*h , , 1,2,..., Step(ii): Initialization Set =y() Step(iii): Recursive approximation Compute successive approximations or the solution using Euler formula =+hF(,), , 1,2,..., -1 Step(iv): Numerical solution The numerical approximation of the solution at 𝑥 = is given by Output: Solution
F(x,y)=
dy
dx
y()=
x
0
y
0
x∈[,]
x
0
x
n
h=-
x
n
x
0
n
n=-
x
n
x
0
h
x
1
x
0
i=0
n
y
0
x
0
y
i+1
y
i
x
i
y
i
i=0
n
𝑥
i+1
y(x)≈y()=
|
x=
x
i+1
x
i+1
y
i+1
y(x)
Below is the Wolfram programming code for Euler method
Below is the Wolfram programming code for Euler method
In[]:=
ClearAll[eulermethodTable]eulermethodTable[f_,start_,end_,y0_,h_]:=Module[{n,x,y},(*Numberofsteps*)n=Round[(end-start)/h];(*x-values*)x=Table[start+ih,{i,0,n}];(*Initializey-valuesarray*)y=ConstantArray[0,n+1];y[[1]]=y0;(*Euleriteration*)Table[y[[i+1]]=y[[i]]+hf[x[[i]],y[[i]]],{i,1,n}];Transpose[{x,y}]]
Runge-Kutta Forth Order
Runge-Kutta Forth Order
The Euler method is derived by truncating the Taylor series expansion after the first derivative term. In contrast, Runge–Kutta methods retain higher-order information from the Taylor expansion, thereby improving numerical accuracy. Consequently, Runge–Kutta methods generally yield more accurate solutions than the Euler method. In this work, attention is given to the classical fourth-order Runge–Kutta method (RK4).For the first-order ODE given in Equation (1) =F(x,y)The fourth-order-Runge-Kutta also known as classical Runge-Kutta or RK4 for the differential equation (1) is given as:
dy
dx
y
i+1
y
i
h
6
k
1
k
2
k
3
k
4
k
1
x
i
y
i
k
2
x
i
h
2
y
i
k
1
2
k
3
x
i
h
2
y
i
k
2
2
k
4
x
i
y
i
k
3
Algorithm: Classical Runge-Kutta method
Algorithm: Classical Runge-Kutta method
Input: Differential equation, with initial condition ; Domain, ; Step size, ; Number of subintervals, . Step(i): Discretization of the domain Define a uniform grid of points =+i*h , , 1,2,..., Step(ii): Initialization Set =y() Step(iii): Recursive approximation For , 1, 2, ..., -1 compute =F(,)=F+,+=F+,+=F(+h,+) and update =+(+2*+2*+), Step(iv): Numerical solution The numerical approximation of the solution at 𝑥 = is given by Output: Solution
F(x,y)=
dy
dx
y()=
x
0
y
0
x∈[,]
x
0
x
n
h=-
x
n
x
0
n
n=-
x
n
x
0
h
x
1
x
0
i=0
n
y
0
x
0
i=0
n
k
1
x
i
y
i
k
2
x
i
h
2
y
i
k
1
2
k
3
x
i
h
2
y
i
k
2
2
k
4
x
i
y
i
k
3
y
i+1
y
i
h
6
k
1
k
2
k
3
k
4
𝑥
i+1
y(x)≈y()=
|
x=
x
i+1
x
i+1
y
i+1
y(x)
Below is the Wolfram programming code for classical Runge-Kutta:
Below is the Wolfram programming code for classical Runge-Kutta:
Formulation of problems
Formulation of problems
Analytical solution
Analytical solution
Implementation of Euler method
Implementation of Euler method
Percentage Error for the Euler method in Exercise 1
Percentage Error for the Euler method in Exercise 1
Solution for the classical Runge-Kutta:
Solution for the classical Runge-Kutta:
Comparison of Euler method, Classical Runge-Kutta, and Analytical solution
Comparison of Euler method, Classical Runge-Kutta, and Analytical solution
Nonlinear first order ordinary differential equations
Exercise 2 : Nonlinear ODE Analysis
Objective: Evaluate the performance of Euler and RK4 methods on a “stiff” nonlinear equation where derivatives change rapidly.
1. Define the equation
2. Define Step size and x - values
Exact solution
Exact solution
3. Exact Solution (NDSolve)
Euler method for exercise 2.
Euler method for exercise 2.
4. Euler Method
Classical Runge - Kutta method for exercise 2.
Classical Runge - Kutta method for exercise 2.
5. RK4 Method
Calculate the errors
Calculate the errors
6. Calculate Absolute Errors
Percentage Error for the Euler method in Exercise 2
Percentage Error for the Euler method in Exercise 2
Construct and display of a comparison table
Construct and display of a comparison table
7. Construct and Display Table
Plot Euler vs. RK4 vs. NDSolveValue
Plot Euler vs. RK4 vs. NDSolveValue
Result and findings
Result and findings
Exercise 2. shows that for nonlinear equations, the Euler method is highly unreliable, producing an error of approximately 39.2 % by the end of the interval . In contrast, the RK4 method maintains a high level of accuracy, staying within 0.2 % of the exact solution .
Exercise 3 : Real - Life Application (Terminal Velocity)
Objective: Apply numerical methods to a physical system: a falling object experiencing atmospheric drag.
1. Define the Physics Equation
Parameters
Parameters
2. Set Parameters
Exact Solution
Exact Solution
3. Exact Solution
Numerical Methods
Numerical Methods
4. Numerical Methods
Create Comparison Table
Create Comparison Table
Plot
Plot
Conclusion
Conclusion
Future work
Future work
While this project focused on first-order equations, there is still more to explore in computational modeling. Future work could look at how these methods handle higher order ODEs. Similarly, it could consider systems of equations, where multiple variables interact simultaneously.
Acknowledgement
Acknowledgement
I sincerely appreciate Wolfram U for providing free access to courses on Introduction to Calculus, Introduction to Differential Equations, and Introduction to Partial Differential Equations, among others. This project was undertaken as a direct outcome of these courses and serves to demonstrate the knowledge and computational skills I have acquired through this learning experience.
References
References
[1] C. Henry Edwards and David E. Penney, Differential Equations & Linear Algebra, 3rd Edition, 2010.
[2] S. Chakraverty, N. R. Mahato, P. Karunakar, and T. D. Rao, Advanced Numerical and Semi-Analytical Methods for Differential Equations. John Wiley & Sons, Inc., 2019.
[3] NASA Glenn Research Center, “Motion of Free Falling Object,” Beginner’s Guide to Aeronautics. Available at: https://www1.grc.nasa.gov/beginners-guide-to-aeronautics/motion-of-free-falling-object/
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
A computational analysis of Euler and classical Runge-Kutta accuracy and stability
by Lateefat Olanike Aselebe
Wolfram Community, STAFF PICKS, February 2, 2026
https://community.wolfram.com/groups/-/m/t/3632288
by Lateefat Olanike Aselebe
Wolfram Community, STAFF PICKS, February 2, 2026
https://community.wolfram.com/groups/-/m/t/3632288