Physics education for the 21st century, solving analytically unsolvable differential equations

The computer revolution in physics

The impact of the computer revolution on physics and engineering can be summarized as follows:

• “Classical physics, in its many branches, is largely formulated in terms of differential equations. Whether we are dealing with mechanics, electromagnetism, or thermodynamics, the behavior of physical systems is described by differential equations.” (1)
• “The vast majority of differential equations that arise in practice from physical, chemical, biological, or engineering problems do not have solutions that can be expressed in terms of elementary functions.” (2)
• “Computational calculus has transformed the way engineers and scientists analyze physical processes. Outside the academic environment, it is now standard practice to use numerical methods to solve problems that were traditionally too complex for analytical approaches.” (3)

The mid-20th century witnessed a silent revolution: simulation based on computational calculus replaced analytic calculus as the primary method for analyzing physical systems in science and engineering outside the classroom. In addition:

• “Euler’s method is one of the simplest and most intuitive approaches to solving differential equations numerically. Its basic concept is trivially easy to grasp., making it an excellent starting point for understanding numerical methods for ODEs.” (4)

The computer revolution marks a fundamental shift in the mathematics of solving differential equations. The power of classical physics increased to the point where now almost all classical systems can be fully analyzed (their performance predicted), with the notable exception of turbulence.

The computer revolution also made it possible to analyze massive engineered systems composed of many interacting subsystems, each with its own differential equation or computational model—something that was inconceivable in the analytic era.

Analytic calculus remains indispensable for proving general theorems—such as the elliptical nature of Keplerian orbits, sinusoidal form of simple harmonic motion, the optimality of the brachistochrone —while the analysis of specific physical systems has shifted almost entirely to computational methods.

The computer revolution is missing in physics education

The computer revolution removed the two limitations of analytic calculus in physics education:

• The differential equation models of most real systems are analytically unsolvable.
• Solving differential equations is an advanced topic in analytic calculus, so instruction in real classical physics cannot begin until the junior year.

The computer revolution completely erases these limitations: computational calculus can compute solutions to all differential equation models of real physical systems, and the basics of computational calculus can be taught in the first week of the first course in physics—in high school or the university—allowing instruction in classical physics to begin immediately (4, 5).

Yet in the current undergraduate curriculum, numerical simulation of differential equation models, if it appears at all, is typically limited to a cursory introduction in an upper-division elective in computational physics. Less frequently, it appears in courses in computational fluid dynamics or computational electrodynamics, which are usually offered at the graduate level.

Academia took note of the computer revolution in physics in Using Computers in Teaching Physics (1989), subtitled, “Computers can revolutionize not only the way we teach physics, but also what physics we teach”. Quoting:

“The computer has revolutionized the way we do physics, but surprisingly, it has not significantly altered the way we teach physics. Talks and papers on teaching with computers fill the meetings and journals of the American Association of Physics Teachers, and workshops on the topic abound, yet the real impact of computers in the classroom is slight.” (6)

Remarkably, nearly 40 years later, little has changed. As noted in Preparing Students for 21st Century Careers (2017):

“Despite the emergence of new disciplinary sub-areas such as computational physics, biophysics, and materials physics, the undergraduate physics curriculum has changed little over the last 50 years.” (7).

And, in Computing in Physics Education (2024):

“In the USA the undergraduate physics curriculum— that is, the standard set of coursework and activities in an accredited physics major—looks much the same as it did 60 years ago.” (8)

Today, simulation is recognized as a useful teaching aid in physics and engineering education. What is far less appreciated is that simulation is the primary method for analyzing differential equation models of physical systems.

Unsolvable systems are invisible in physics education

While the word unsolvable rarely, if ever, appears in physics textbooks or physics education literature, unsolvable systems are occasionally alluded to.

From Classical Mechanics:

“The three Euler equations determine the motion of a spinning body as seen in a frame fixed in the body. In general, they are difficult to use because the components Γ₁, Γ₂, Γ₃ of the applied torque as seen in the rotating body frame are complicated (and unknown) functions of time. In fact, the main use of Euler’s equations is in the case that the applied torque is zero.” (9)

Euler’s equations are the basis for the analysis of 3-D rigid body motion; when the applied torques are non-zero the equations are analytically unsolvable.

From Classical Dynamics:

“The actual motion of a rocket attempting to leave Earth’s gravitational field is quite complicated. For analytical purposes, we begin by making several assumptions. The rocket will have only vertical motion, with no horizontal component. We neglect air resistance and assume the acceleration of gravity is constant. We also assume that the burn rate of the fuel is constant. All these factors that are neglected can reasonably be included with a numerical analysis by computer”. (1)

From Heat and Mass Transfer:

“Analytical solution methods are limited to highly simplified problems in simple geometries. The geometry must be such that its entire surface can be described mathematically in a coordinate system by setting the variables equal to constants.”

“Even minor complications in geometry can make an analytical solution impossible. Even in simple geometries, heat transfer problems cannot be solved analytically if the thermal conditions are not sufficiently simple.”

“Numerical solution of engineering problems has now become the norm rather than the exception even when analytical solutions are available.” (10)

From Mechanics of Materials:

“Apart from a few cases of simple beam geometry and loading, this equation is difficult to solve, because it represents a nonlinear second-order differential equation. Fortunately, it can be modified to….. ꝺ²v/ꝺx² = M/El.” (11)

From Fundamentals of Fluid Mechanics:

“Unfortunately, because of the general complexity of the Navier–Stokes equations (they are nonlinear, second-order, partial differential equations), they are not amenable to exact mathematical solutions except in a few instances.” (12)

From Numerical solution of initial boundary problems involving Maxwell’s equations in isotropic media:

“Solutions to the time-dependent Maxwell’s equations in general form are unknown except for a few special cases.” (13)

How are unsolvable systems covered in undergraduate physics and engineering education? The answer is: they aren’t, they aren’t mentioned. The limitations of analytic calculus have been internalized in the current curriculum to the extent that unsolvable systems are invisible.

Incorporating simulation into physics and engineering education

The mathematical basis for simulation is computational calculus. Unlike analytic calculus, which requires several semesters of dedicated mathematics coursework, computational calculus is conceptually simple and intuitively transparent. The basics of Euler’s method—the foundation of computational calculus—can be taught to high-school science students with no prior exposure to calculus in a single one-hour lecture (4, 5).

The Coming Revolution in Physics Education (5) describes a high school course in which students use computational calculus to analyze real systems, e.g., the Apollo trajectory, modeled with ordinary differential equations (ODEs). This paper extends that approach to simulate real systems modeled by partial differential equations (PDEs), e.g., an electromagnetic wave hitting a barrier with an opening.

Together, the two papers contain analyses of eight prototypical systems spanning classical physics that can form the basis for a two-semester introductory course in classical physics. Figure 1 shows the results of these analyses.

The simplicity and transparency of computational calculus

Simulation erases the limitations of analytic calculus, and solving differential equation models of physical systems, which is usually difficult or impossible analytically, becomes easy computationally.

Take, for example, Newton’s falling apple. Newton discovered the laws of gravity and motion, and derived the governing differential equation for a falling object, a = G·m_E/r². This equation is derived in high school physics. However, the equation is very difficult to solve analytically, and solving it is beyond the scope of not only of high school physics, but also current university physics (14). A representative analytic solution is given in Wikipedia Free fall (15), in the form of an infinite series:

y ( t ) = ∑ n = 1 ∞ x n n ! lim q → 0 ( d n − 1 d r n − 1 [ r n ( 7 2 ( arcsin ( r ) − r − r 2 ) ) − 2 3 n ] )

Evaluating this expression yields.

y ( t ) = y 0 ( x − 1 5 x 2 − 3 175 x 3 − 23 7875 x 4 − 1894 3931875 x 5 − 3293 21896875 x 6 − 2418092 62077640625 x 7 − ⋯ ) , x = [ 3 2 ( π 2 − t 2 u y 0 3 ) ] 2 3

The Coming Revolution in Physics Education (5) shows how students in a high school physics class used Euler’s method to simulate a wide variety of systems modeled by ordinary differential equations, including falling objects, orbits, and rocket trajectories.

The differential equation model for Newton’s falling apple has two state variables, the velocity of the apple v and the position of the apple r, i.e., its distance from the center of the earth. The differential equations in the model are:

r ’ ( t ) = v ( t )

v ’ ( t ) = G · m E / r ( t ) 2

The computational loop for a falling apple follows immediately from the model, it is:for i=1:N % Loop for Euler’s computational equations. t(i+1) = t(i) + dt; % Advance time. A = −GmE/ r(i)^2; % Compute acceleration. r(i+1) = r(i) + v(i)dt; % Eq. (1). v(i+1) = v(i) + A*dt; % Eq. (2). end

The computational solution really is this simple and intuitively transparent. High school students can use trigonometry to resolve acceleration into its x and y components, and calculate 2-D orbits and rocket trajectories, e.g., the Apollo trajectory above (5). This is how it’s done in science and engineering today, using Runge–Kutta methods, which are refinements of Euler’s method (3).

Systems modeled by partial differential equations

Euler’s method for ODEs is easily extended to compute solutions to PDEs, and is known as the finite difference method (FDM); the translation from the partial differential equations in the system model to arithmetic computational equations is simple, 1-to-1, intuitively clear, and by rote.

Suppose T(t,x) is a function for the temperature T in a wire at time t and position x. Given a partial differential equation for T of the form.

∂ T ( t , x ) / ∂ t = c · ∂ 2 T ( t , x ) / ∂ x 2 ,

the FDM uses this approximation for the 2nd order partial derivative ∂²T(t,x)/∂x²

∂ 2 T ( t , x ) ∂ x 2 ≈ ∂ T ( t , x + Δ x / 2 ) / ∂ x − ∂ T ( t , x − Δ x / 2 ) / ∂ x Δ x

≈ T ( t , x + Δ x ) − T ( t , x ) Δ x − T ( t , x ) − T ( t , x − Δ x ) Δ x Δ x

= T ( t , x + Δ x ) − 2 · T ( t , x ) − T ( t , x − Δ x ) Δ x 2

Using Euler’s formula for the time derivative, the computational equation for T is:

T ( t + Δ t , x ) ≈ T ( t , x ) + c ⋅ ∂ T ( t , x ) ∂ t ⋅ Δ t = T ( t , x ) + c ⋅ ∂ 2 T ( t , x ) ∂ x 2 ⋅ Δ t ≈ T ( t , x ) + c ⋅ T ( t , x + Δ x ) − 2 · T ( t , x ) + T ( t , x − Δ x ) Δ x 2 ⋅ Δ t

The FDM is the only math used in this paper to solve differential equations.

System analyses

Since the time of Newton, the basic paradigm for the analysis of physical systems has been:

• State the laws of physics governing the system using differential equations.
• Derive a differential equation model of the system from the laws in step 1.
• Analyze the model, with the goal of predicting the performance of the system.

Differential equation models are concise descriptions of the physics of physical systems, but they do not enable us to directly predict the behavior of that system. Euler’s method and the FDM translate the differential equations 1-to-1 to computational equations which are also concise and do enable us to precisely predict the behavior of the system.

The mathematics required for the system analyses below consists of only the FDM. The system models are derived in Appendices I and II using the analytic definitions of first and second order derivatives and contour integrals; no additional analytic calculus is required, not even the product rule or the chain rule. The contour integrals are evaluated only for infinitesimal rectangular contours, so they can be approximated by their corresponding finite difference expressions. (Tensors and a little more calculus appear in Appendices IV and V).

Heat conductance

Computing 1-D time-dependent heat conductance

The partial differential equation model for 1-D time-dependent heat conductance (Appendix I):

∂ T ( t , x ) / ∂ t = ( k / c ρ ) ⋅ ( ∂ T 2 ( t , x ) ) / ( ∂ x 2 ) where k is conductivity, c is heat capacity, ρ is density

The time domain 0-T is divided into N evenly spaced intervals 0 = t₁, … t_{*N* + 1} = T, with Δt = T/N, and the spatial domain 0-X is divided into a uniform grid of points, with 0 = x₁, … x_{*Nx* + 1} = X, with ∆x = X/N_X.

Using the FDM,

T ( t n + 1 , x i ) = T ( t n , x i ) + ( k c ρ ) ⋅ ( ∂ 2 T ( t n , x i ) ∂ x 2 ) ⋅ Δ t

Making the FDM substitution for the 2nd order partial derivative gives the computational equation, where the state variables are the temperatures at the grid points x_i, 1 ≤ I ≤ N + 1:

T ( t n + 1 , x i ) = T ( t n , x i ) + ( k c ρ ) ⋅ ( T ( t n , x i + 1 ) − 2 ⋅ T ( t n , x i ) + ( T ( t n , x i − 1 ) ) ⋅ Δ t

Figure 2 shows 1-D heat conductance with 11 grid points, Δx = 1, with starting temperatures T = [0 5 5 5 5–5 –5 –5 –5 –5 2.5] and the endpoints held to their initial fixed values

Computing 2-D time-dependent heat conductance

The partial differential equation model for 2-dimensional time-dependent heat conductance (Appendix I):

∂ T ( t , x , y ) ∂ t = ( k c ρ ) ⋅ ( ∂ 2 T ( t , x , y ) ∂ x 2 + ∂ 2 T ( t , x , y ) ∂ y 2 )

To analyze heat conductance in a 2-D rectangular plate the time domain 0-T is divided into N evenly spaced intervals 0 = t₁, … t_{*N* + 1} = T, with Δt = T/N, and the spatial domain 0-X, 0-Y is divided into a uniform grid of points as shown in the diagram, with 0 = x₁, … x_{*NX* + 1} = X, with ∆x = X/N_X, and 0 = y₁, … y_{*Ny* + 1} = Y, with ∆y = Y/N_Y.

The state variables are the temperature values T(x_i, y_j) at the grid points.

The computational equation for T is:

T n + 1 ( x i , y j ) = T n ( x i , y j ) + ( k c ρ ) ⋅ ( ∂ 2 T n ( x i , y j ) ∂ x 2 + ∂ 2 T n ( x i , y j ) ∂ y 2 ) ⋅ Δ t

following the standard convention of writing Tⁿ(x_i,y_j) for T(t_n,x_i,y_j).

Making the FDM substitutions for the partial derivatives:

T n + 1 ( x i , y j ) = T n ( x i , y j ) + ( k c ρ ) ⋅ ( ( ( T n ( x i + 1 , y j ) − 2 T n ( x i , y j ) + T n ( x i − 1 , y j ) ) / Δ x 2 + ( T n ( x i , y j + 1 ) − 2 T n ( x i , y j ) + T n ( x i , y j − 1 ) ) / Δ y 2 ) ⋅ Δ t )

With ∆x = ∆y this simplifies to

T n + 1 ( x i , y j ) = T n ( x i , y j ) + ( k c ρ ) ⋅ ( T n ( x i + 1 , y j ) + T n ( x i , y j + 1 ) − 4 T n ( x i , y j ) + T n ( x i − 1 , y j ) + T n ( x i , y j − 1 ) Δ x 2 ) ⋅ Δ t

The above equation written in index form is

T n + 1 ( x i , y j ) = T n ( x i , y j ) + ( k c ρ ) ⋅ ( T n ( x i + 1 , y j ) + T n ( x i , y j + 1 ) − 4 T n ( x i , y j ) + T n ( x i − 1 , y j ) + T n ( x i , y j − 1 ) Δ x 2 ) ⋅ Δ t

Given the values for Tⁿ(i,j), i.e., the temperature values for the grid at time t_n, the computational equations can be used to calculate T^{*n* + 1}(i,j) for each of the interior points in the grid.

The grid has four boundaries, the left, Tⁿ(1,j), right, Tⁿ(N_X + 1,j) bottom, Tⁿ(i,1), top, Tⁿ(i, N_Y + 1).

There are three common ways of specifying boundary values, as follows:

• Dirichlet boundary condition—the boundary value is constant, i.e., does not change with time, e.g., T^{*n* + 1}(i,1) = C_i.
• Neumann boundary condition—the value of a partial derivative is specified at the boundary, e.g., (T^{*n* + 1}(i,2) − T^{*n* + 1}(i,1))/Δy = S_i, so T^{*n* + 1}(i,1) = T^{*n* + 1}(i,2)∙S_i∙Δy.
• Robin boundary condition—a linear combination of Dirichlet and Neumann boundary conditions.

It’s necessary to specify initial conditions, i.e., temperatures, for all of the grid points at time t₁ = 0. In the example shown in Figure 3 Dirichlet boundary conditions keep the bottom boundary at 0 °K, the upper boundary at 5 °K, and the interior section at 10 °K, Neumann boundary conditions ∂T/∂x = 0 are specified for the left and right borders.

Simulations

1-D heat eq—https://www.youtube.com/watch?v=VOgzw9Cbvfk.

2-D heat eq—https://www.youtube.com/watch?v=G0wt_qIa6II.

Electrodynamics

Computing 1-D EM waves

Maxwell’s 3rd and 4th equations in component form are the model for electromagnetic radiation. In 2-D, with I = 0 they are, where μ is magnetic permeability and ε is electric permittivity (Appendix II):

∂ E y ∂ x − ∂ E x ∂ y = − ∂ B z ∂ t

∂ E z ∂ y − ∂ E y ∂ z = − ∂ B x ∂ t

∂ E x ∂ z − ∂ E z ∂ x = − ∂ B y ∂ t

∂ B y ∂ x − ∂ B x ∂ y = μ ε ∂ E z ∂ t

∂ B z ∂ y − ∂ B y ∂ z = μ ε ∂ E x ∂ t

∂ B x ∂ z − ∂ B z ∂ x = μ ε ∂ E y ∂ t

In 1-D, we start with an E field that is in the z direction and varies with x, thus E_z(x) represents our E field and E_x(x) = E_y(x) = 0. From the 1st set of equations above:

∂ E y ∂ x − ∂ E x ∂ y = − ∂ B z ∂ t = 0

∂ E z ∂ y − ∂ E y ∂ z = − ∂ B x ∂ t = 0

∂ E z ∂ x = ∂ B y ∂ t

Thus B_z and B_x are constant and assumed 0. From the 2nd set of equations, we have

∂ B y ∂ x = μ ε ∂ E z ∂ t

Rewriting the model in terms of H and E (with H = B/μ):

∂ H y ∂ t = 1 μ ∂ E z ∂ x

∂ E z ∂ t = 1 ε ∂ H y ∂ x

To compute solutions to the equations we will make the usual FDM substitutions with one modification, instead of using the Euler approximation for the 1st order time derivative, we will use the more accurate ‘leapfrog’ central difference approximation (in anticipation of the Yee algorithm):

∂ H y ∂ t ≈ H y n + 1 ( i ) − H y n − 1 ( i ) 2 ⋅ Δ t

So we have

H y n + 1 ( i ) = H y n − 1 ( i ) + 2 ⋅ Δ t μ ⋅ E z n ( i + 1 ) − E z n ( i − 1 ) 2 ⋅ Δ x

and

E z n + 1 ( i ) = E z n − 1 ( i ) + 2 ⋅ Δ t ε ⋅ H y n ( i + 1 ) − H y n ( i − 1 ) 2 ⋅ Δ x

To scale ∆t and ∆x we set the Courant number SC = c ∙ ∆t / ∆x = 1, where c is the speed of light. Thus ∆t is commensurate with the time it takes for the wave to move from one grid point to the next, i.e. ∆x, and the equations become:

H y n + 1 ( i ) = H y n − 1 ( i ) + 1 c ⋅ μ ⋅ ( E z n ( i + 1 , j ) − E z n ( i − 1 ) )

and similarly,

E z n + 1 ( i ) = E z n − 1 ( i ) + 1 c ⋅ ε ⋅ ( H y n ( i + 1 ) − H y n ( i − 1 ) )

Figure 4 shows a 1-D electromagnetic wave originating at x = 0 and traveling in free space and encountering a dielectric with relative permittivity ε_R = 5 at x = 300 Δx. when it hits the dielectric surface the wave is partially reflected (note the E field is inverted) and partially transmitted.

Computing 2-D EM waves

We start with an E field that is in the z direction and varies with x and y, thus E_z(x,y) represents our E field, and E_x(x,y) = E_y(x,y) = 0. From the first set of equations, we have (with I = 0 and H=B/):

− μ ∂ H x ∂ t = ∂ E z ∂ y

μ ∂ H y ∂ t = ∂ E z ∂ x

And from the 2nd set

ε ∂ E z ∂ t = ∂ H y ∂ x − ∂ H x ∂ y

Making the FDM substitutions yields the computational equations, with SC =1 and ∆x = ∆y:

H x n + 1 ( i , j ) = H x n − 1 ( i , j ) + 1 c ⋅ μ ⋅ ( E z n ( i + 1 , j ) − E z n ( i − 1 , j ) )

H y n + 1 ( i , j ) = H y n − 1 ( i , j ) + 1 c ⋅ μ ⋅ ( E z n ( i + 1 , j ) − E z n ( i − 1 , j ) )

E z n + 1 ( i , j ) = E z n − 1 ( i , j ) + 1 c ⋅ ε ⋅ ( H y n ( i + 1 , j ) − H y n ( i − 1 , j ) − H x n ( i + 1 , j ) − H x n ( i − 1 , j ) )

A perfect electric conductor (PEC) has zero resistance and infinite conductivity. An electric field cannot exist in a perfect electric conductor. A perfect electric conductor on a boundary is modeled by a Dirichlet boundary E_z(t,x,y) = 0. The magnetic field is reflected and the corresponding Neumann boundary condition is ∂H_y(t,x,y)/∂x = 0.

Figure 5 shows a 2-D wave before and after hitting a barrier with an opening. A Ricker wavelet is generated on the left boundary, at 40 Δx in it hits a PEC barrier that has an opening in the center, the Neumann boundary conditions for the top and bottom boundaries are ꝺEz/ꝺy = ꝺHx/ꝺy = ꝺHy/ꝺy = 0.

Simulations

1-D EM wave—https://www.youtube.com/watch?v=tIOo8vqwdfA

2-D EM wave—https://www.youtube.com\watch?v=PaCSUhlUTps

The Yee algorithm

The state-of-the-art computational method for electromagnetic waves is the Yee algorithm (13). The Yee algorithm can be viewed (and programmed) as a restricted version of the leapfrog FDM.

Figure 6 shows a 1-D FDM grid with spacing Δx and Δt. If we consider the subset of points (n, i) with both n and i even or both n and i odd, we have a staggered Yee grid with spacing 2Δx and 2Δt. If we propagate the wave on that subset using the Yee algorithm the calculations are identical to the calculations for propagating the wave using the leapfrog FDM. That is, Yee algorithm produces the same values as the leapfrog FDM for these points. The 1-D Yee algorithm calculates ½*½ as many values as the 1-D leapfrog FDM, and is hence 4 times faster, with identical accuracy but less resolution.

Similarly, if we start with a 2-D FDM grid with spacing Δx, Δx, Δt and calculate the solution on the corresponding staggered Yee grid with spacing 2Δx, 2Δx, and 2Δt, i.e. E_zⁿ(i,j) for n, i, and j all odd, H_xⁿ(i,j) for n even, i odd and j even, H_y ⁿ(i,j) for n even, i even, and j odd, the Yee solution is identical to the leapfrog FDM solution for these points. The 2-D Yee algorithm calculates ½½½ as many values as the 2-D leapfrog FDM and is hence 8 times faster, with identical accuracy but less resolution.

Summary

The computer revolutionized physics and engineering in the mid-20th century by making it possible to solve virtually all differential equation models of physical systems. Yet the anticipated revolution in physics and engineering education has not happened. The curriculum remains constrained by the invisible limitations of analytic calculus, even though computational calculus has long been the primary method for analyzing physical systems in science and industry.

This paper, together with The Coming Revolution in Physics Education, demonstrates that realistic physical systems can be studied from the very first physics course. In the proposed introductory course:

• Differential equations are used to model physical systems. This is the defining starting point of the course.
• Euler’s method is used to simulate systems modeled by ordinary differential equations, and the finite difference method is used for systems modeled by partial differential equations.
• Students implement their own simulations, building both modeling skills and computational competence.

This is possible in an introductory course because computational calculus has two extraordinary features:

• Conceptually, it is trivially easy and is intuitively transparent; the basic yet powerful Euler’s method can be taught to high school science students in a single one-hour lecture.
• Practically, it can solve all differential equation models of real physical systems.

Simulation, when built by students from first principles, enables them to visualize observable system behavior and to graph the evolution of invisible system and auxiliary variables. The connections between the governing laws of physics, the differential equation model, and the computational are explicit and transparent. The computational equations demystify the differential equation model because they show exactly how each system variable is updated on a millisecond-by-millisecond basis. In addition, students can easily explore “what if” questions by varying system or environmental parameters, running the simulation, and receiving immediate feedback. This approach represents a new dimension in physics and engineering education.

Together, the two papers contain complete computational analyses of prototypical systems spanning classical physics. These analyses can form the basis for a two-semester freshman course. Students starting their sophomore year will have a strong foundation in the analysis of differential equation models of physical systems. This will have a profound impact on the rest of the STEM curriculum. A new curriculum incorporating computational calculus will give students a sense of accomplishment and mastery that the current curriculum does not provide, and align physics and engineering education with physics and engineering practice in the 21st century.

Generative AI statement

The author(s) declared that Generative AI was used in the creation of this manuscript. Various editing details were discussed with generative AI.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2026.1884678/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Nenad Manojlović, University of Algarve, Portugal

Reviewed by: Imtiaz Ahmad, University of Swabi, Pakistan

Hikmawati Hikmawati, University of Mataram, Indonesia