Numerical Differentiation

There are many ways to evaluate the differentiation of a function
numerically, one of which is three-point method. Due to its method which
require three point to evaluate an derivative, first and last point required
a different form of formula which given below.

$$\begin{array}{}\text{(1)}& f^{\prime }(x_0)={\displaystyle \frac{1}{2h}}[-3f(x_0)+4f(x_0+h)-f(x_0+2h)]+{\displaystyle \frac{h^2}{3}}{f}^{(3)}({\zeta }_1)\end{array}$$

$$\begin{array}{}\text{(2)}& f^{\prime }(x_0)={\displaystyle \frac{1}{2h}}[f(x_0-2h)-4f(x_0-h)+3f(x_0)]+{\displaystyle \frac{h^2}{3}}{f}^{(3)}({\zeta }_2)\end{array}$$

$$\begin{array}{}\text{(3)}& f^{\prime }(x_0)={\displaystyle \frac{1}{2h}}[f(x_0+h)+-f(x_0-h)]+{\displaystyle \frac{h^2}{6}}{f}^{(3)}({\zeta }_3)\end{array}$$
By intuition, the information on both end is less than the one in the middle,
so you can see from the error function that the error is two times larger
at boundary than the one in the middle.
Using the error function which is give by the last term from (1) to (3), we can
know complete the table which compute the derivative of two different function.

Example 1

$$f(x)=e^{2x}$$

Example 2

$$f(x)=xlnx$$

Code

program threepoint
!========================================================
!   Purpose: Find derivative
!
!   Methods: Three point method
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!     10/30/20       MorrisH        Original Code
!========================================================
    implicit none
    double precision, external :: f, g
    double precision :: x(4) = (/ 8.1d0, 8.3d0, 8.5d0, 8.7d0 /)
    double precision :: x1(4), h=.2d0, err(4), ebd(4)
    integer :: i

    ! calculate derivative
    do i = 1, 4
        if (i == 1) then
            x1(i) = ( -3.d0 * f(x(i)) + 4.d0 * f(x(i+1)) - &
                      f(x(i+2)) ) / (2.d0 * h)
        else if (i == 4) then
            x1(i) = ( f(x(i-2)) - 4.d0 * f(x(i-1)) + &
                      3.d0 * f(x(i)) ) / (2.d0 * h)
        else
            x1(i) = (f(x(i+1)) - f(x(i-1))) / (2.d0 * h)
        end if
    end do

    ! calculate error
    do i = 1, 4
        err(i) = abs(x1(i) - g(x(i))) / g(x(i))
    end do

    ! calculate error bound
    do i = 1, 4
        if ((i==1) .or. (i==4)) then
            ebd(i) = abs(h**2 / 3.d0 * 8.d0 * (-1.d0 / x(i)**2))
        else
            ebd(i) = abs(h**2 / 6.d0 * 8.d0 * (-1.d0 / x(i)**2))
        end if
    end do

    ! print the table
    write(6, '(5A15)') 'x', 'f(x)', "f'(x)", 'Error', 'Error Bound'
    do i = 1, 4
        write(6, '(5F15.5)') x(i), f(x(i)), x1(i), err(i), ebd(i)
    end do
end program threepoint

function f(x)
    implicit none
    double precision :: f, x
    ! f = dexp(2.d0 * x)
    f = x * log(x)
end function

function g(x)
    implicit none
    double precision :: g, x
    ! g = 2.d0 * dexp(2.d0 * x)
    g = log(x) + 1.d0
end function

Example 3: Circuit problem

A circuit with impressed voltage

$\epsilon (t)$ and inductance

$L$,
can be express as following equation.

$$\begin{array}{}\text{(4)}& \epsilon =L{\displaystyle \frac{di}{dt}}+Ri\end{array}$$
Again, using three-point method to approximate the impressed volage given by (4).
With

$L=0.98$ and

$R=0.142$ the result table is given below.

Code

program rl
!========================================================
!   Purpose: Solve RL circuit
!
!   Methods: Three point method
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!     10/30/20       MorrisH        Original Code
!========================================================
    implicit none
    double precision :: t(5) = (/ 1.d0, 1.01d0, 1.02d0, 1.03d0, 1.04d0 /)
    double precision :: x(5) = (/ 3.1d0, 3.12d0, 3.14d0, 3.18d0, 3.24d0 /)
    double precision :: x1(5), v(5), h=.01d0
    integer :: i

    ! calculate derivative
    do i = 1, 5
        if (i == 1) then
            x1(i) = ( -3.d0 * x(i) + 4.d0 * x(i+1) - &
                      x(i+2) ) / (2.d0 * h)
        else if (i == 5) then
            x1(i) = ( x(i-2) - 4.d0 * x(i-1) + &
                      3.d0 * x(i) ) / (2.d0 * h)
        else
            x1(i) = (x(i+1) - x(i-1)) / (2.d0 * h)
        end if
    end do

    ! calculate voltage
    do i = 1, 5
        v(i) = .98d0 * x1(i) + .142d0 * x(i)
    end do

    ! print the table
    write(6, '(4A15)') 't', 'i', "di/dt", 'voltage'
    do i = 1, 5
        write(6, '(4F15.5)') t(i), x(i), x1(i), v(i)
    end do
end program rl

Example 4: Predator-Prey model

Predator-Prey model is govern by following two differential equation
where

$x_1(x)$ is the number of preys and

$x_2(x)$ is the number of
predators.

$$\begin{array}{rl}{\displaystyle \frac{d{x}_1(t)}{dt}}& =k_1{x}_1-k_2{x}_1{x}_2\\ {\displaystyle \frac{d{x}_2(t)}{dt}}& =k_3{x}_1{x}_2-k_4{x}_2\end{array}$$
Solving these equation using RK4 method accompany with initial value value given below
gives the following result.

$$\begin{array}{rl}& x_1(0)=1000\\ & x_2(0)=500\\ & k_1=3\\ & k_2=0.002\\ & k_3=0.0006\\ & k_4=0.5\end{array}$$
Due to the number of predators initally is outnumber by preys,
number of prey raise up. But it didn't last long, soon it follow
by the raising of predator which cause the number of prey decrease
dracmatically and almost distinguise. And this cycle goes on forever
both predator and prey never settle to a single number.

Code

program sysrk4
!========================================================
!   purpose: solve system of ode
!
!   methods: rk4
!       
!       date        programer       description of change
!       ====        =========       =====================
!     10/29/20       morrish        original code
!========================================================
    implicit none
    double precision, external :: f1, f2
    double precision :: k(4), w(4, 2), h, t
    double precision, dimension(:), allocatable :: x1, x2
    integer :: i, n

    ! coef
    k = (/ 3.d0, 2.d-3, 6.d-4, .5d0 /)
    ! time step
    h = 1.d-3
    n = nint(20/h)
    allocate(x1(n))
    allocate(x2(n))
    ! inital value
    x1(1) = 1000.d0
    x2(1) = 500.d0

    do i = 1, n
        w(1, 1) = h * f1(x1(i), x2(i), k)
        w(1, 2) = h * f2(x1(i), x2(i), k)

        w(2, 1) = h * f1(x1(i)+0.5d0*w(1, 1), &
                         x2(i)+0.5d0*w(1, 2), k)
        w(2, 2) = h * f2(x1(i)+0.5d0*w(1, 1), &
                         x2(i)+0.5d0*w(1, 2), k)

        w(3, 1) = h * f1(x1(i)+0.5d0*w(2, 1), &
                         x2(i)+0.5d0*w(2, 2), k)
        w(3, 2) = h * f2(x1(i)+0.5d0*w(2, 1), &
                         x2(i)+0.5d0*w(2, 2), k)

        w(4, 1) = h * f1(x1(i)+w(3, 1), x2(i)+w(3, 2), k)
        w(4, 2) = h * f2(x1(i)+w(3, 1), x2(i)+w(3, 2), k)

        x1(i+1) = x1(i) + (w(1, 1) + 2.d0 * w(2, 1) + &
                           2.d0 * w(3, 1) + w(4, 1)) / 6.d0
        x2(i+1) = x2(i) + (w(1, 2) + 2.d0 * w(2, 2) + &
                           2.d0 * w(3, 2) + w(4, 2)) / 6.d0
    end do

    open(66, file='out.dat', status='unknown')
    do i = 1, n
        write(66, *) dble(i-1)*h, x1(i), x2(i)
    end do
end program sysrk4


function f1(x1, x2, k)
    implicit none
    double precision :: f1, x1, x2
    double precision :: k(4)
    f1 = k(1) * x1 - k(2) * x1 * x2 
end function


function f2(x1, x2, k)
    implicit none
    double precision :: f2, x1, x2
    double precision :: k(4)
    f2 = k(3) * x1 * x2 - k(4) * x2
end function

Example 5: Projectile motion

Consider a baseball was being hit by 2020 World Series MVP,
Corey Seager, we can write down the equation of motion of the
baseball by

$$\mathbf{F}=m\mathbf{g}-{\displaystyle \frac{1}{2}}{c}_w\rho A{v}^2\hat{v}$$
where

$c_w$ is the dimensionless drag coefficient,

$\rho$ is
the air density,

$v$ is the velocity of the ball, and

$A$ is
the cross-sectional area of the ball.

Without air resistance

We first ignore resistance, we can see from below figure
that the inital speed for it to clear the fence is roughly
around 25.5 m/s.

As you can see, this system should always satisfied energy
conservation, but it can loose accaracy when preform numerical
calculation. So far so good, using RK4 energy were conserved
which also can be a indicator to the converges of the parameter
that we use in RK4 method(eg. time step).

Following is the result of that calculation, it show that the
initial speed for the baseball to clear the fence with 35 degrees
angle is around 25.432810 m/s.

Code

program proj
!========================================================
!   purpose: solve system of ode
!
!   methods: rk4
!       
!       date        programer       description of change
!       ====        =========       =====================
!     10/29/20       morrish        original code
!========================================================
    implicit none
    double precision, parameter :: pi = 3.1415926535898d0
    double precision, external :: f1, f2, f3, f4
    double precision :: w(4, 4), h, theta, v0
    double precision, dimension(:), allocatable :: x, y, vx, vy
    double precision, dimension(:), allocatable :: t
    integer :: i, j, n
    character(len=15) :: str

    ! time step
    h = 1.d-2
    n = nint(100/h)
    allocate(x(n))
    allocate(y(n))
    allocate(vx(n))
    allocate(vy(n))
    allocate(t(n))

    ! inital value
    theta = 35.d0 * pi / 180.d0
    v0 = 25.0

    speed: do 
        t = 0.d0
        x = 0.d0
        y = 0.d0
        vx = 0.d0
        vy = 0.d0
        y(1) = .7d0
        vx(1) = v0 * dcos(theta)
        vy(1) = v0 * dsin(theta)
        do i = 1, n
            w(1, 1) = h * f1(x(i), y(i), vx(i), vy(i))
            w(1, 2) = h * f2(x(i), y(i), vx(i), vy(i))
            w(1, 3) = h * f3(x(i), y(i), vx(i), vy(i))
            w(1, 4) = h * f4(x(i), y(i), vx(i), vy(i))

            w(2, 1) = h * f1(x(i)+0.5d0*w(1, 1), &
                             y(i)+0.5d0*w(1, 2), &
                             vx(i)+0.5d0*w(1, 3), &
                             vy(i)+0.5d0*w(1, 4))
            w(2, 2) = h * f2(x(i)+0.5d0*w(1, 1), &
                             y(i)+0.5d0*w(1, 2), &
                             vx(i)+0.5d0*w(1, 3), &
                             vy(i)+0.5d0*w(1, 4))
            w(2, 3) = h * f3(x(i)+0.5d0*w(1, 1), &
                             y(i)+0.5d0*w(1, 2), &
                             vx(i)+0.5d0*w(1, 3), &
                             vy(i)+0.5d0*w(1, 4))
            w(2, 4) = h * f4(x(i)+0.5d0*w(1, 1), &
                             y(i)+0.5d0*w(1, 2), &
                             vx(i)+0.5d0*w(1, 3), &
                             vy(i)+0.5d0*w(1, 4))

            w(3, 1) = h * f1(x(i)+0.5d0*w(2, 1), &
                             y(i)+0.5d0*w(2, 2), &
                             vx(i)+0.5d0*w(2, 3), &
                             vy(i)+0.5d0*w(2, 4))
            w(3, 2) = h * f2(x(i)+0.5d0*w(2, 1), &
                             y(i)+0.5d0*w(2, 2), &
                             vx(i)+0.5d0*w(2, 3), &
                             vy(i)+0.5d0*w(2, 4))
            w(3, 3) = h * f3(x(i)+0.5d0*w(2, 1), &
                             y(i)+0.5d0*w(2, 2), &
                             vx(i)+0.5d0*w(2, 3), &
                             vy(i)+0.5d0*w(2, 4))
            w(3, 4) = h * f4(x(i)+0.5d0*w(2, 1), &
                             y(i)+0.5d0*w(2, 2), &
                             vx(i)+0.5d0*w(2, 3), &
                             vy(i)+0.5d0*w(2, 4))

            w(4, 1) = h * f1(x(i)+w(3, 1), y(i)+w(3, 2), &
                             vx(i)+w(3, 3), vy(i)+w(3, 4))
            w(4, 2) = h * f2(x(i)+w(3, 1), y(i)+w(3, 2), &
                             vx(i)+w(3, 3), vy(i)+w(3, 4))
            w(4, 3) = h * f3(x(i)+w(3, 1), y(i)+w(3, 2), &
                             vx(i)+w(3, 3), vy(i)+w(3, 4))
            w(4, 4) = h * f4(x(i)+w(3, 1), y(i)+w(3, 2), &
                             vx(i)+w(3, 3), vy(i)+w(3, 4))

            x(i+1) = x(i) + (w(1, 1) + 2.d0 * w(2, 1) + &
                             2.d0 * w(3, 1) + w(4, 1)) / 6.d0
            y(i+1) = y(i) + (w(1, 2) + 2.d0 * w(2, 2) + &
                             2.d0 * w(3, 2) + w(4, 2)) / 6.d0
            vx(i+1) = vx(i) + (w(1, 3) + 2.d0 * w(2, 3) + &
                               2.d0 * w(3, 3) + w(4, 3)) / 6.d0
            vy(i+1) = vy(i) + (w(1, 4) + 2.d0 * w(2, 4) + &
                               2.d0 * w(3, 4) + w(4, 4)) / 6.d0
            t(i+1) = dble(i)*h
            ! hit the groud
            if (y(i+1) <= 0.d0) then
                exit
            end if
        end do
        ! check if it pass the fence
        do j = 1, i+1 
            if ((x(j) >= 60.d0).and.(y(j) >= 2.d0)) then
                exit speed
            end if
        end do
        ! write to a file
        v0 = v0 + 0.00001d0
        ! write(str, '(I6.6)') j
        ! open(66, file='data/a/'//trim(str)//'.dat', &
        !      status='unknown')
        ! do i = 1, n
        !     write(66, *) t(i), x(i), y(i), vx(i), vy(i)
        !     if (y(i) <= 0.d0) then
        !         exit
        !     end if
        ! end do
        ! close(66)
    end do speed
    write(6, *) 'With inital angle of 35 degrees, the velocity when &
                 the ball pass the fence.'
    write(6, '(5A15)') 'v0', 'x(final)', 'y(final)', 'vx(final)', 'vy(final)'
    write(6, '(5f15.6)') v0, x(j), y(j), vx(j), vy(j)

end program proj


function f1(x, y, vx, vy)
    implicit none
    double precision :: f1, x, y, vx, vy
    f1 = vx
end function

function f2(x, y, vx, vy)
    implicit none
    double precision :: f2, x, y, vx, vy
    f2 = vy
end function

function f3(x, y, vx, vy)
    implicit none
    double precision :: f3, x, y, vx, vy
    f3 = 0.d0
end function

function f4(x, y, vx, vy)
    implicit none
    double precision :: f4, x, y, vx, vy
    f4 = -9.8d0
end function
\end{minted}

With air reistance

Next we take a step further, considering air resistance
by changing differenct initial velocity I found out that
in order to clear the fence we need at least 35.57340 m/s.

With this speed in mind we try to find the best possible
angle that can easily clear the fence. We can see from
the trajectories, the best angle is roughly around 40
degrees.

Below figure verified our guess, when the angle is
40.550 we can vertically clear the fence by 1.112 m.

program proj
!========================================================
!   purpose: solve system of ode
!
!   methods: rk4
!       
!       date        programer       description of change
!       ====        =========       =====================
!     10/29/20       morrish        original code
!========================================================
    implicit none
    double precision, parameter :: pi = 3.1415926535898d0
    double precision, external :: f1, f2, f3, f4
    double precision :: w(4, 4), h, theta, v0
    double precision, dimension(:), allocatable :: x, y, vx, vy
    double precision, dimension(:), allocatable :: t
    integer :: i, j, n
    character(len=15) :: str

    ! time step
    h = 1.d-3
    n = nint(100/h)
    allocate(x(n))
    allocate(y(n))
    allocate(vx(n))
    allocate(vy(n))
    allocate(t(n))

    ! inital value
    theta = 30.d0 * pi / 180.d0
    v0 = 35.573740d0

    speed: do 
        t = 0.d0
        x = 0.d0
        y = 0.d0
        vx = 0.d0
        vy = 0.d0
        y(1) = .7d0
        vx(1) = v0 * dcos(theta)
        vy(1) = v0 * dsin(theta)
        do i = 1, n
            w(1, 1) = h * f1(x(i), y(i), vx(i), vy(i))
            w(1, 2) = h * f2(x(i), y(i), vx(i), vy(i))
            w(1, 3) = h * f3(x(i), y(i), vx(i), vy(i))
            w(1, 4) = h * f4(x(i), y(i), vx(i), vy(i))

            w(2, 1) = h * f1(x(i)+0.5d0*w(1, 1), &
                             y(i)+0.5d0*w(1, 2), &
                             vx(i)+0.5d0*w(1, 3), &
                             vy(i)+0.5d0*w(1, 4))
            w(2, 2) = h * f2(x(i)+0.5d0*w(1, 1), &
                             y(i)+0.5d0*w(1, 2), &
                             vx(i)+0.5d0*w(1, 3), &
                             vy(i)+0.5d0*w(1, 4))
            w(2, 3) = h * f3(x(i)+0.5d0*w(1, 1), &
                             y(i)+0.5d0*w(1, 2), &
                             vx(i)+0.5d0*w(1, 3), &
                             vy(i)+0.5d0*w(1, 4))
            w(2, 4) = h * f4(x(i)+0.5d0*w(1, 1), &
                             y(i)+0.5d0*w(1, 2), &
                             vx(i)+0.5d0*w(1, 3), &
                             vy(i)+0.5d0*w(1, 4))

            w(3, 1) = h * f1(x(i)+0.5d0*w(2, 1), &
                             y(i)+0.5d0*w(2, 2), &
                             vx(i)+0.5d0*w(2, 3), &
                             vy(i)+0.5d0*w(2, 4))
            w(3, 2) = h * f2(x(i)+0.5d0*w(2, 1), &
                             y(i)+0.5d0*w(2, 2), &
                             vx(i)+0.5d0*w(2, 3), &
                             vy(i)+0.5d0*w(2, 4))
            w(3, 3) = h * f3(x(i)+0.5d0*w(2, 1), &
                             y(i)+0.5d0*w(2, 2), &
                             vx(i)+0.5d0*w(2, 3), &
                             vy(i)+0.5d0*w(2, 4))
            w(3, 4) = h * f4(x(i)+0.5d0*w(2, 1), &
                             y(i)+0.5d0*w(2, 2), &
                             vx(i)+0.5d0*w(2, 3), &
                             vy(i)+0.5d0*w(2, 4))

            w(4, 1) = h * f1(x(i)+w(3, 1), y(i)+w(3, 2), &
                             vx(i)+w(3, 3), vy(i)+w(3, 4))
            w(4, 2) = h * f2(x(i)+w(3, 1), y(i)+w(3, 2), &
                             vx(i)+w(3, 3), vy(i)+w(3, 4))
            w(4, 3) = h * f3(x(i)+w(3, 1), y(i)+w(3, 2), &
                             vx(i)+w(3, 3), vy(i)+w(3, 4))
            w(4, 4) = h * f4(x(i)+w(3, 1), y(i)+w(3, 2), &
                             vx(i)+w(3, 3), vy(i)+w(3, 4))

            x(i+1) = x(i) + (w(1, 1) + 2.d0 * w(2, 1) + &
                             2.d0 * w(3, 1) + w(4, 1)) / 6.d0
            y(i+1) = y(i) + (w(1, 2) + 2.d0 * w(2, 2) + &
                             2.d0 * w(3, 2) + w(4, 2)) / 6.d0
            vx(i+1) = vx(i) + (w(1, 3) + 2.d0 * w(2, 3) + &
                               2.d0 * w(3, 3) + w(4, 3)) / 6.d0
            vy(i+1) = vy(i) + (w(1, 4) + 2.d0 * w(2, 4) + &
                               2.d0 * w(3, 4) + w(4, 4)) / 6.d0
            t(i+1) = dble(i)*h
            ! hit the groud
            if (y(i+1) <= 0.d0) then
                exit
            end if
        end do
        ! check if it pass the fence
        open(66, file='data/angle.dat', status='unknown')
        do j = 1, i+1 
            if (x(j) >= 60.d0) then
                write(66, *) theta / pi * 180.d0, y(j)
                exit
            end if
        end do
        theta = (theta / pi * 180.d0 + .001d0) * pi / 180.d0
        if (theta / pi * 180.d0 >= 50.d0) then 
            exit speed
        end if

        ! write to a file
        ! write(str, '(I6.6)') j
        ! open(66, file='data/c/'//trim(str)//'.dat', &
        !      status='unknown')
        ! do i = 1, n
        !     write(66, *) t(i), x(i), y(i), vx(i), vy(i)
        !     if (y(i) <= 0.d0) then
        !         exit
        !     end if
        ! end do
        ! close(66)
    end do speed
    ! write(6, *) 'With inital speed of 35.57374 m/s, the angle when &
    !              the ball pass the fence.'
    ! write(6, '(5A15)') 'angle', 'x(final)', 'y(final)', 'vx(final)', 'vy(final)'
    ! write(6, '(5f15.6)') theta/pi*180.d0, x(j), y(j), vx(j), vy(j)

end program proj


function f1(x, y, vx, vy)
    implicit none
    double precision :: f1, x, y, vx, vy
    f1 = vx
end function

function f2(x, y, vx, vy)
    implicit none
    double precision :: f2, x, y, vx, vy
    f2 = vy
end function

function f3(x, y, vx, vy)
    implicit none
    double precision :: f3, x, y, vx, vy
    f3 = -0.5d0 * .5d0 * 1.335d0 * 7.8539d-3 * vx * &
         dsqrt(vx**2 + vy**2) / .2d0
end function

function f4(x, y, vx, vy)
    implicit none
    double precision :: f4, x, y, vx, vy
    f4 = -9.8d0 - 0.5d0 * .5d0 * 1.335d0 * 7.8539d-3 * vy * &
         dsqrt(vx**2 + vy**2) / .2d0

end function