Here are some excerpts taken from notes I wrote on the topic when I was a student:
The Fundamental Theorem of Calculus (FTOC) provides the vital connecting link between the two tools of calculus: differentiation and integration. Use of this link offers a way to demonstrate that the methods for approximating definite integrals and the numerical methods for approximating solutions to first order ordinary initial value problems are mathematically equivalent. To illustrate this, let's begin with the so-called derivative form of the FTOC:
The existence of an indefinite integral of $\displaystyle f(x)$, i.e., of a function $\displaystyle F(x)$ whose derivative is $\displaystyle f(x)$ is given by:
(1) $\displaystyle \frac{d}{dx}\int_a^x f(u)\,du=f(x)$ or $\displaystyle F(x)+C=\int_a^x f(u)\,du$, C a constant.
It should be remarked that the variable u in the symbol $\displaystyle f(u)\,du$ given in (1) is only a "dummy variable." For the variable of integration, whatever it may be called, is "integrated out" and disappears. Having established this, let's go on to define:
$\displaystyle F(x_1)+C=\int_a^{x_1} f(u)\,du$ and $\displaystyle F(x_2)+C=\int_a^{x_2} f(u)\,du$ where $\displaystyle x_1<x_2$
Then we have:
$\displaystyle \left(F(x_1)+C \right)-\left(F(x_2)+C \right)=\left(\int_a^{x_1} f(u)\,du \right)-\left(\int_a^{x_2} f(u)\,du \right)$
$\displaystyle F(x_1)-F(x_2)=\left(\int_a^{x_1} f(u)\,du \right)-\left(\int_a^{x_1} f(u)\,du+\int_{x_1}^{x_2} f(u)\,du \right)$
(2) $\displaystyle \int_{x_1}^{x_2} f(u)\,du=F(x_2)-F(x_1)$
This is the anti-derivative form of the FTOC, and from this we may also conclude:
$\displaystyle \int_a^b f(u)\,du=-\int_b^a f(u)\,du$
We have tacitly assumed that $\displaystyle f(x)$ is continuous and differentiable on $\displaystyle [a,b]$.
While the existence of indefinite and definite integrals of continuous functions is established, the technique of finding them may be far from simple. In fact in many cases the integrals of elementary functions cannot be expressed in terms of elementary functions themselves. For example, consider the simple functions:
$\displaystyle f(x)=\sqrt{x^3+1}$ and $\displaystyle f(x)=e^{x^2}$
While there are substitutions and techniques to transform many into forms found in a finite table of integrals, the evaluation of a definite integral may prove exceedingly difficult if not impossible. It therefore becomes useful to to develop methods for approximating definite integrals as they find a wide application in physical problems. Besides the use in plane areas, the definite integral is used for volumes, arc lengths, surface areas, rectilinear motion, center of mass, moments of inertia, electrostatic and gravitational potentials, liquid pressure, biology, economics, etc.
Approximate Integration
Suppose we have some approximation or numeric scheme A such that:
(4) $\displaystyle A\approx\int_a^b f(x,y)\,dx$, where $\displaystyle \frac{dy}{dx}=f(x,y)$
Since the slope of y is assumed to be continuous on the interval $\displaystyle [a,b]$, given some initial value for y, we may suppose that an explicit relationship between x and y can be found within the rectangle $\displaystyle (a,y(a))-(b,y(b))$ by the uniqueness and existence theorem.
Now, also suppose we have divided the interval $\displaystyle [a,b]$ into n sub-intervals of equal width where:
$\displaystyle a\le x_k<x_{k+1}\le b$ and $\displaystyle x_0=a,\,x_n=b,\,0\le k\le n-1$.
Since by the additive property we have:
$\displaystyle \int_a^b f(x,y)\,dx=\sum_{k=0}^{n-1}\left[\int_{x_k}^{x_{k+1}} f(x,y)\,dx \right]$ we may define:
(5) $\displaystyle A=\sum_{k=0}^{n-1}A_k$ where $\displaystyle A_k\approx\int_{x_k}^{x_{k+1}} f(x,y)\,dx$
If the approximation converges, then we must have:
$\displaystyle \lim_{n\to\infty}\left[\sum_{k=0}^{n-1}A_k \right]=\int_a^b f(x,y)\,dx$
By defining $\displaystyle y_k\equiv y(x_k)$ and applying (2) we have:
$\displaystyle A_n\approx y_{n+1}-y_n$
Solving for $\displaystyle y_{n+1}$ gives rise to the numeric scheme:
(6) $\displaystyle y_{n+1}\approx y_n+A_n$
Thus, the approximating summation for the definite integral has become an approximation to the first order IVP:
(6a) $\displaystyle \frac{dy}{dx}=f(x,y)$, $\displaystyle y(x_0)=y_0$
We now have a means of approximating the solution $\displaystyle y(x)$ at $\displaystyle x=x_n$ where $\displaystyle a\le x_n\le b$.
Riemann Sums and the Approximation Method of Euler
In my opinion, one of the most straight-forward approximation methods available for definite integrals are the Riemann sums with regular partitions, which approximates the definite integral with a series of rectangles of equal width and whose heights are the function's value at $\displaystyle x_n$. But, as usual, the tradeoff for the simplicity of the method is that it does not converge very rapidly.
First, approximate the following definite integral using a Riemann sum of one partition:
(10) $\displaystyle \int_{x_n}^{x_{n+1}} f(x,y)\,dx\approx\Delta x\cdot f(x,y)$, where $\displaystyle \Delta x=x_{n+1}-x_n$
which leads to the numeric scheme:
(11) $\displaystyle y_{n+1}\approx y_n+\Delta x\cdot f(x,y)$, $\displaystyle y_0=y(x_0)$
where $\displaystyle A_n=\Delta x\cdot f(x,y)$ which is base $\displaystyle \Delta x$ times height $\displaystyle f(x,y)$. This is an example of an explicit scheme, i.e., it can be solved for $\displaystyle y_{n+1}$.
Approximation Method of Euler (tangent line method)
When we use the direction field method to sketch a particular solution to an IVP, we try to visualize the intermediate directions between the isoclines we have drawn. If we follow a finite number of these directions, the sketch becomes a polygonal curve or chain of line segments. This polygonal curve is, visually speaking, an approximation to the solution. We can construct values $\displaystyle y_n$ that approximate the solution values $\displaystyle y(x_n)$ as follows:
One method we may use to demonstrate the derivation of Euler's method is through the use of the differential to obtain a linear approximation (the tangent line). Another method would be to use the point-slope formula or Taylor formula of order 1. At the point $\displaystyle (x_n,y_n)$, the slope of the solution is given by:
$\displaystyle \frac{dy}{dx}=f(x_n,y_n)$
Recall that the definition of the differential of the dependent variable is $\displaystyle \Delta y\approx\Delta x\frac{dy}{dx}$
Using $\displaystyle \Delta y=y_{n+1}-y_{n}$ this yields the recursive formula:
(12) $\displaystyle y_{n+1}\approx y_n+\Delta x\cdot f(x,y)$
Thus, the solution at $\displaystyle x=x_n$ may be approximated by (12). This is equivalent to (11), showing that the Riemann sums and approximation method of Euler are equivalent.
In a similar fashion, the trapezoidal method and the improved Euler's method are analogs, as are the mid-point rule and the second order Runge-Kutta method, and Simpson's Rule and the fourth order Runge-Kutta method, but I will leave these for later. (Happy)