Suppose we want to find 

x₄
òf(x)dx
x₀

Many different numerical methods exist to approximate such an integral. Here we use a 4-point Newton-Cotes formula.

We can approximate f(x) in the interval [x₀,x₄] by a polynomial, then analytically integrate the polynomial term by term. The higher the degree of the polynomial, the greater the accuracy of the resulting approximation to the integral. Here we use the fourth degree polynomial f(x)=a x⁴+b x³+c x²+d x+e, and use the method of Undetermined Coefficients to find a, b, c, d and e. If we then analytically integrate this polynomial term by term, the result is an approximation to the integral of f(x) from x₀ to x₄, and it is:

x₄
òf(x)dx»2h/45[7f(x₀)+32f(x₁)+12f(x₂)+32f(x₃)+7f(x₄)]
x₀

where h=(x₄-x₀)/4 and the selected points are equally spaced.  If we have n subintervals, with n being a multiple of 4, we have the following extended formula,

x_n
òf(x)dx»2h/45[7f(x₀)+32f(x₁)+12f(x₂)+32f(x₃)+14f(x₄)+...+7f(x_n)]
x₀

in which h=(x_n-x₀)/n:

The integration is performed once and the result stored.  It is then performed again with twice as many subintervals and the result compared with the previous one.  The process is continued until two conditions are met:

a) The result differs from the previous one by less than the specified error.

b) The number of subintervals is at least the number specified by the user.

MULTIPLE INTEGRATION

Suppose what we want is a Multiple Integral. OK, then suppose f(x) above is in turn an integral, for example, over the variable y, with limits that may be functions of x.  We would then have a double integral

 d     b(x)
ò  ò g(x,y) dy dx
 c     a(x)

If we group it like this:

 d            b(x)
ò [  ò g(x,y) dy ]dx
 c            a(x)

we see that the part inside the brackets is something like f(x).

We would first perform the integration over y. For such integration over y we would in turn use the method previously described. Again, y itself may be an integral over another variable, say z, and so on.

Each integration is performed once and the result stored.  It is then performed again with twice as many subintervals and the result compared with the previous one.  The process is continued until two conditions are met:

a) The result differs from the previous one by less than the specified error.

b) The number of subintervals is at least the number specified by the user.

DIFFERENTIATION

Some central difference formulas to find derivates:

First derivative:

               f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)
f ’(x) = ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                       12h

Second derivative:

               -f(x-2h)+16f(x-h)-30f(x)+16f(x+h)-f(x+2h)
f ’’(x) = ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                              12h²

Higher derivatives:

    1           n                         n!
¾¾¾     S   (-1)ⁿ  ¾¾¾¾¾¾   f( x+ (n-2k)h )
(2h)ⁿ     k=0               (n-k)!  k!

In all cases h is a small quantity.

PARTIAL DIFFERENTIATION

Here we use the convention

          ¶F                   ¶²F                   ¶²F
F_x = ¾¾  , F_xx = ¾¾¾  , F_xy = ¾¾¾
          ¶x                   ¶x²                  ¶x¶y

and similarly for all other combinations.

Suppose we have a function of 4 variables F(x,y,z,t). If h is a small quantity, then the first and second order partial derivatives of F(x,y,z,t) can be approximated by the following expressions. Here we assume that the function and all of its 1st and 2nd order partial derivatives are continuous at the point we’re interested in.

            F(x-2h,y,z,t)-8*F(x-h,y,z,t)+8*F(x+h,y,z,t)-F(x+2h,y,z,t)
F_x =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                            12h

 
            F(x,y-2h,z,t)-8*F(x,y-h,z,t)+8*F(x,y+h,z,t)-F(x,y+2h,z,t)
F_y =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                            12h

           F(x,y,z-2h,t)-8*F(x,y,z-h,t)+8*F(x,y,z+h,t)-F(x,y,z+2h,t)
F_z =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                            12h

           F(x,y,z,t-2h)-8*F(x,y,z,t-h)+8*F(x,y,z,t+h)-F(x,y,z,t+2h)
F_t =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                           12h

              -F(x-2h,y,z,t)+16*F(x-h,y,z,t)-30*F(x,y,z,t)+16*F(x+h,y,z,t)-F(x+2h,y,z,t)
F_xx =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                               12h²

              -F(x,y-2h,z,t)+16*F(x,y-h,z,t)-30*F(x,y,z,t)+16*F(x,y+h,z,t)-F(x,y+2h,z,t)
F_yy =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                               12h²

             -F(x,y,z-2h,t)+16*F(x,y,z-h,t)-30*F(x,y,z,t)+16*F(x,y,z+h,t)-F(x,y,z+2h,t)
F_zz =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                              12h²

             -F(x,y,z,t-2h)+16*F(x,y,z,t-h)-30*F(x,y,z,t)+16*F(x,y,z,t+h)-F(x,y,z,t+2h)
Ftt =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                              12h²

               F_x(x,y-2h,z,t)-8*F_x(x,y-h,z,t)+8*F_x(x,y+h,z,t)-F_x(x,y+2h,z,t)
F_xy =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                  12h

               F_x(x,y,z-2h,t)-8*F_x(x,y,z-h,t)+8*F_x(x,y,z+h,t)-F_x(x,y,z+2h,t)
F_xz =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                 12h

               F_x(x,y,z,t-2h)-8*F_x(x,y,z,t-h)+8*F_x(x,y,z,t+h)-F_x(x,y,z,t+2h)
F_xt =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                 12h

               F_y(x,y,z-2h,t)-8*F_y(x,y,z-h,t)+8*F_y(x,y,z+h,t)-F_y(x,y,z+2h,t)
F_yz =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                 12h

               F_y(x,y,z,t-2h)-8*F_y(x,y,z,t-h)+8*F_y(x,y,z,t+h)-F_y(x,y,z,t+2h)
F_yt =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                 12h

               F_z(x,y,z,t-2h)-8*F_z(x,y,z,t-h)+8*F_z(x,y,z,t+h)-F_z(x,y,z,t+2h)
F_zt =  ¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾¾
                                                                 12h

From these, it’s easy to generalize to any number of independent variables.