Numerical Differentiation and Integration

Differentiation

The derivative of the function $f$ at $x_0$ is

\[f'(x_0) = lim_{h \to 0}\frac{f(x_0+h) - f(x_0)}{h}\]

Three and Five Point formula.

when $x_1 = x_0+h$ and $x_2 = x_0 + 2h$ for some $h\ne0$ gives

Three-Point Formulas

\[f'(x_0) = \frac{1}{h}[-\frac{3}{2}f(x_0) + 2f(x_1) - \frac{1}{2}f(x_2)] + \frac{h^2}{3}f^{(3)}(\xi_0)\]

Three-Point Endpoint Formula

\[f'(x_0) = \frac{1}{2h}[-3f(x_0)+4f(x_0+h) - f(x_0+2h)] + \frac{h^2}{3}f^{(3)}(\xi_0)\]

Three-Point Midpoint formula

\[f'(x_0) = \frac{1}{2h}[f(x_0+h) - f(x_0-h)] - \frac{h^2}{6}f^{(3)}(\xi_0)\]

using NumericalAnalysis: NCalculus
NCalculus.ThreePoint(x->sin(x), 0, 0.1)

1.0033216789612571

Five-Point MidPoint formula

\[f'(x_0) = \frac{1}{12h}[f(x_0 - 2h) - 8f(x_0 -h)+8f(x_0+h)-f(x_0 + 2h)]+\frac{h^4}{30}f^{(5)}(\xi_0)\]

Five-Point EndPoint formula

\[f'(x_0) = \frac{1}{12h}[-25f(x_0)+48f(x_0 +h)-36f(x_0+2h)+16f(x_0+3h)-3f(x_0+4h)]+\frac{h^4}{5}f^{(5)}(\xi_0)\]

NCalculus.FivePoint(x->sin(x), 0, 0.1)

0.9999803084008577

Integration

Trapezoidal Rule

\[\int_a^b f(x)dx = \frac{h}{2}[f(a)+f(b)] - \frac{h^3}{12}f''(\xi)\]

NCalculus.Trapezoidal(x->x^2, 0, 2)

4.0

Simpson's Rule

\[\int_a^b f(x)dx = \frac{h}{3}[f(a) + f(a+h) + f(b)] - \frac{h^5}{90}f^{(4)}(\xi)\]

NCalculus.Simpson(x->x^2, 0, 2)

2.6666666666666665

Newton-Cotes Formulas

will calculate n=1,2,3,4. formula rule. if you input n in $[1,4]$, will get a value.

NCalculus.Newton_cotes(x->sin(x),0,π/4)

4-element Array{Float64,1}:
 0.30055886494217315
 0.29798754218726264
 0.2928586591925902
 0.2928692281360844

NCalculus.Newton_cotes(x->sin(x), 0, π/4, n=4)

0.2928692281360844

NCalculus.Newton_cotes(x->sin(x),0,π/4,n=3, method="Open")

0.29291070254917145

Romberg Integration

\[R_{k,j} = R_{k,j-1} + \frac{1}{4^{j-1} - 1}(R_{k,j-1} - R_{k-1, j-1}), k = j,j+1,...\]

NCalculus.Romberg(x->sin(x), 0, π/4, 5, seq_tab=true)

5×5 Array{Float64,2}:
 0.27768   0.0       0.0       0.0       0.0
 0.28912   0.292933  0.0       0.0       0.0
 0.291952  0.292896  0.292893  0.0       0.0
 0.292658  0.292893  0.292893  0.292893  0.0
 0.292834  0.292893  0.292893  0.292893  0.292893

NCalculus.Romberg(x->sin(x), 0, π/4, 5)

0.29289321881345187

Gaussian Quadrature

Suppose we want to determine $c_1, c_2, x_1$ and $x_2$ so

\[\int_{-1}^{1}f(x)dx=c_1f(x_1) + c_2f(x_2)\]

we use a Ploynomial $f(x) = a_0 + a_1x+a_2x^2+a_3x^3$

we get: a solution $c_1 = c_2=1, x_1 = -\frac{\sqrt{3}}{3}, x_2=\frac{\sqrt{3}}{3}$

So:

\[\int_{-1}^{1}f(x)dx \approx f(\frac{-\sqrt{3}}{3}) + f(\frac{\sqrt{3}}{3})\]

NCalculus.Gaussian_Quad(x->sin(x),n=5, a=0, b=π/4)

0.29289321879993296

NCalculus.Gaussian_Quad(x->x^2)

0.6666666666485936

Multiple Integrals

In this section, Consider the double integral

\[\iint_R f(x,y)dA\]

where $R = {(x,y) | a \leq x \leq b, c \leq y \leq d}$. for some constants a, b, c, d, is a rectangular in the plane.

we apply the Composite Trapezoidal rule to calculate:

Simpson Double Integral

\[\int_a^b\int_{ydown(x)}^{yup(x)}f(x, y)dydx\]

f(x, y) = log(x + 2y)
y_up(x) = 1.5
y_down(x) = 1.0
NCalculus.SimpsonDoubleIntegral(f, (1.4, 2.0), y_up, y_down)

0.4295544977526275

Gaussian Double Integral

NCalculus.GaussianDoubleIntegral(f, (1.4, 2.0), y_up, y_down)

0.42955452750549294

Method

NumericalAnalysis.NCalculus — Module

NCalculus

Numerical Differentiation and Integration Module.

source

NumericalAnalysis.NCalculus.FivePoint — Method

FivePoint(f::Function, x₀::Real, h::Real; method::String="Endpoint")

input a function f, x₀ and h.Args is method (default: "Endpoint", if you set method is't "Endpoint", it will return MidPoint ans.)

source

NumericalAnalysis.NCalculus.GaussianDoubleIntegral — Method

 GaussianDoubleIntegral(f::Function, interval_x::Tuple, y_up::Function, y_down::Function, m::Int=5, n::Int=5)

$f = f(x, y), x ∈ interval_x, y = yup(x) - ydown(x)$, Args m and n is Separation distance.

source

NumericalAnalysis.NCalculus.Gaussian_Quad — Method

Gaussian_Quad(f::Function; n::Int=3, a::Real=-1, b::Real=1)

Gaussian Quadrature must set n<6,default calculate int f(x)dx in [-1,1]. you can set interval [a, b] to change default interval.

source

NumericalAnalysis.NCalculus.MutipleIntegralSimple — Method

MutipleIntegralSimple(f::Function, interval_x::Tuple, interval_y::Tuple)

It will apply the Composite Trapezoidal rule to approximate the integral.

source

NumericalAnalysis.NCalculus.Newton_cotes — Method

Newton_Cotes(f::Function, a::Real, b::Real; n::Int=0, method::String="Midpoint")

input a function f, and $x in [a, b]$calculateint_a^b f(x)dx. argsn∈[1,4]`

source

NumericalAnalysis.NCalculus.Romberg — Method

Romberg(f::Function, a::Real, b::Real, n::Int; seq_tab::Bool=false)

input a function f, and $x in [a, b]$, set n.if Args seq_tab=true, will return a table, else return a value.

source

NumericalAnalysis.NCalculus.Simpson — Method

Simpson(f::Function, a::Real, b::Real)

input a function f, and $x in [a, b]$ calculate ∫f(x)dx.

source

NumericalAnalysis.NCalculus.SimpsonDoubleIntegral — Method

SimpsonDoubleIntegral(f::Function, interval_x::Tuple, y_up::Function, y_down::Function; n::Int =6, m::Int=6)

$f = f(x, y), x ∈ interval_x, y = yup(x) - ydown(x)$, Args m and n is Separation distance.

source

NumericalAnalysis.NCalculus.ThreePoint — Method

ThreePoint(f::Function, x₀::Real, h::Real; method::String="Endpoint")

input a function f, x₀ and h.Args is method (default: "Endpoint", if you set method is't "Endpoint", it will return MidPoint ans.)

source

NumericalAnalysis.NCalculus.Trapezoidal — Method

Trapezoidal(f::function, a::Real, b::Real)

input a function f, and $x in [a, b]$ calculate ∫f(x)dx.

source