Why is it called numerical integration when we numerically solve differential equations? This has been bugging me literally for years.
When numerically simulating a system of differential equations (e.g., with Runge-Kutta or Euler methods), we are using the derivative to estimate the value of the function at the next time step. Why is this called numerical integration or integration rather than simply numerical simulation or function estimation or something?
I have not found this nomenclature discussed, and would love to see the origins. I am probably Googling wrong and missing something obvious. My guess is that from the fundamental theorem, an operation that brings you from $\dot{x}$ to x is by definition integration, so we are technically doing numerical integration?
I suppose it could be that all of the same terms are involved as when you calculate the integral using a Riemann sum (or related techniques). But for the differential equation we are not calculating the area but the value of the function at the next time step, so it doesn't seem like an integral in that sense.
 A: Integration is the general term for the resolution of a differential equation.
You probably know the simple case of antiderivatives,
$$\int f(x)\,dx$$ which in fact solve the ODE $$y'(x)=f(x)$$ via an integral.
The same term is used when you solve, say
$$y'(x)=y(x)+5,$$
giving
$$y(x)=ce^x+5.$$
You integrate the equation. Sometimes, the solution itself is called an integral.
You can integrate by analytical methods, and also by numerical methods.
A: A solver which solves differential equations by integrating variables step by step is called differential equation solver or integrator.
Suppose that it is a forced mass damper spring equation:
Fx=m* d2x+b* dx+k* x , where x is a function of t => x(t), d2x is second derivative of x, dx is first derivative of x with respect to t. k is stiffness coefficient , b is damping coeff. , m is the mass.
We suppose time interval between t[n] and t[n+1] is fixed and it is dt.
create initial values for x(0)=0 and dx(0)=0
You find d2x for next point like this:
d2x[i+1]=(Fx-bdx[i]-kx[i])/m
Then you integrate (accumulate values) Hence the name integrator:
dx[i+1]=dx[i]+d2x[i+1]*dt
x[i+1]=x[i]+dx[i]*dt
then you increase i , i=i+1 and do it all over again ...
This is why it is called integrator...
Here it will be more clearer directly using integrals :
See this section ?
  Fxb+=Fx[i]*dt # first integral of Fx wrt dt 
  Fxb2+=Fxb*dt # second integral of Fx wrt dt**2
  xb+=x[i]*dt # first integral of x wrt dt 
  xb2+=xb*dt # second integral of Fx wrt dt**2

Code:
from matplotlib import pyplot as plt 
from matplotlib import animation, rc
import sympy as sp
import numpy as np
from scipy.integrate import odeint
import time

m=1
b=0.1
k=100
t=sp.symbols('t')
x=sp.Function('x')
Fx=100*sp.sin(3*t)
expr=m*sp.diff(x(t),t,2)+b*sp.diff(x(t),t)+k*x(t)-Fx
x=sp.dsolve(expr,x(t)).rhs
C2=sp.solve(x.subs(t,0))[0]
x=x.subs('C2',C2)
C1=sp.solve(sp.diff(x,t).subs(t,0))[0]
x=x.subs('C1',C1)
print(x)

f=sp.lambdify(t,x)

T=5
dt=0.01
t_analytic=np.linspace(0,T,int(T/dt))
x_analytic=f(t_analytic)

######################################
T=5
dt=0.01
t=np.linspace(0,T,int(T/dt))
m=1
b=0.1
k=100

Fx=100*np.sin(3*t)

x=np.zeros(len(t))

tdump=time.time()
Fxb,Fxb2,xb,xb2=[0,0,0,0]

for i in range(0,len(t)-1):
  Fxb+=Fx[i]*dt # first integral of Fx wrt dt 
  Fxb2+=Fxb*dt # second integral of Fx wrt dt**2
  xb+=x[i]*dt # first integral of x wrt dt 
  xb2+=xb*dt # second integral of Fx wrt dt**2
  x[i+1]=(Fxb2-b*xb-k*xb2)/m

plt.plot(t[2:],x[2:],'black')

print('Integrating solv. Elapsed:',time.time()-tdump)
#plt.plot(t,np.diff(x)/dt,'g--')

err=np.sum(np.abs(x-f(t)))
total=np.sum(np.abs(x_analytic))
print('Integrating solv. error ratio is:',err/total)

def mydiff(x, t):
  m = 1 # Mass
  b = 0.1 # Damping constant
  k = 100 # Stiffness of the spring
  F = 100*np.sin(3*t)

  dx1dt = x[1]
  dx2dt = (F - b*x[1] - k*x[0])/m
  dxdt = [dx1dt, dx2dt]
  return dxdt

T=5
#dt=0.01
t=np.linspace(0,T,int(T/dt))
x_init = [0,0]
# Solve ODE
tdump=time.time()
x = odeint(mydiff, x_init, t)
x1 = x[:,0]
x2 = x[:,1]
print('odeint Elapsed:',time.time()-tdump)
# Plot the Results
plt.plot(t,x1,'r--')

err=np.sum(np.abs(x1-f(t)))
total=np.sum(np.abs(x_analytic))
print('odeint Error ratio is:',err/total)


There are a lot sophisticated version of this itteration. Runge-Kutta , ADAMS method. There implicit and explicit versions of this.
There is an MIT course notes of R.Sureshkumar:
https://web.mit.edu/10.001/Web/Course_Notes/Differential_Equations_Notes/lec24.html
