Minimal working example: NKTg law (varying inertia, RK4, CSV output)

4

Hello @NKTgLaw, I think I see a potential problem in your approach, you have taken up the case of varying inertia but that is precisely the case where we cannot take a = dv/dt = F_{net}/m because

p = mv \\ \frac{dp}{dt} = m\frac{dv}{dt} + v\frac{dm}{dt} \\ \frac{dv}{dt} = \frac{1}{m}\frac{dp}{dt} - \frac{v}{m}\frac{dm}{dt} = \frac{F_{net}}{m} - \frac{v}{m}\frac{dm}{dt}

Or in terms of your notation

\frac{dv}{dt} = \frac{F_{net}}{m} - \frac{NTKg_2}{m^2}

That said, I tried your code too, noticed a few errors with gfortran tweaked it a little and run it for the case of a simple harmonic oscillator, and plotted the following graph using gnuplot

plot
plot640×480 21.1 KB

which is in agreement with theory but if I update the mdot in the source code, it doesn’t change the trajectory of the oscillator at all, which it ideally should.

The modified source code for getting the data points for the above graph is attached below.

modified source code for the datapoints of the above graph 
! nktg_demo.f90
! Minimal working example for the “NKTg Law on Varying Inertia”
! State: x (position), v (velocity), m (mass)
! p = mv; NKTg1 = xp; NKTg2 = (dm/dt)*p
! Dynamics (example):
! dx/dt = v
! dv/dt = (Fx + Fspring + Fdamp + Fthrust)/m
! dm/dt = mdot (negative for fuel burn)
!
! Compile:
! gfortran -O2 -std=f2008 nktg_demo.f90 -o nktg_demo
! Run:
! ./nktg_demo > out.csv

module nktg_types
    implicit none
    private
    integer, parameter, public :: dp = kind(1.0d0)

    type, public :: Params
        real(dp) :: k = 0.0_dp ! spring constant (N/m)
        real(dp) :: c = 0.0_dp ! linear damping (N*s/m)
        real(dp) :: Fx = 0.0_dp ! external constant force (N)
        real(dp) :: thrust = 0.0_dp ! constant thrust (N)
        real(dp) :: mdot = 0.0_dp ! mass rate (kg/s), <0 for burn
    end type Params
contains
end module nktg_types

module nktg_model
    use nktg_types
    implicit none
    private
    public :: rhs, nktg1, nktg2, momentum

contains

    pure subroutine rhs(t, y, par, dydt)
        !! y = [x, v, m]
        real(dp), intent(in) :: t
        real(dp), intent(in) :: y(3)
        type(Params), intent(in) :: par
        real(dp), intent(out) :: dydt(3)

        real(dp) :: x, v, m, p, Fspring, Fdamp, Fnet

        x = y(1);  v = y(2);  m = max(y(3), 1.0e-12_dp)  ! avoid divide-by-zero
        p = m*v

        Fspring = -par%k * x
        Fdamp   = -par%c * v
        Fnet    = par%Fx + Fspring + Fdamp + par%thrust

        dydt(1) = v
        dydt(2) = Fnet / m
        dydt(3) = par%mdot

    end subroutine rhs

    pure real(dp) function momentum(v, m) result(p)
        real(dp), intent(in) :: v, m
        p = m*v
    end function momentum

    pure real(dp) function nktg1(x, v, m) result(val)
        real(dp), intent(in) :: x, v, m
        val = x * momentum(v,m)
    end function nktg1

    pure real(dp) function nktg2(mdot, v, m) result(val)
        real(dp), intent(in) :: mdot, v, m
        val = mdot * momentum(v,m)
    end function nktg2

end module nktg_model

module rk4
    use nktg_types
    implicit none
    private
    public :: step_rk4

contains

    subroutine step_rk4(t, y, h, par, f, yout)
        interface
            pure subroutine f(t, y, par, dydt)
                use nktg_types
                implicit none
                real(kind=kind(1.0d0)), intent(in) :: t
                real(kind=kind(1.0d0)), intent(in) :: y(3)
                type(Params), intent(in) :: par
                real(kind=kind(1.0d0)), intent(out) :: dydt(3)
            end subroutine
        end interface

        real(dp), intent(in) :: t, h
        real(dp), intent(in) :: y(3)
        type(Params), intent(in) :: par
        real(dp), intent(out) :: yout(3)

        real(dp) :: k1(3), k2(3), k3(3), k4(3), yt(3)

        call f(t, y, par, k1)
        yt = y + 0.5_dp*h*k1
        call f(t + 0.5_dp*h, yt, par, k2)
        yt = y + 0.5_dp*h*k2
        call f(t + 0.5_dp*h, yt, par, k3)
        yt = y + h*k3
        call f(t + h, yt, par, k4)

        yout = y + (h/6.0_dp)*(k1 + 2.0_dp*k2 + 2.0_dp*k3 + k4)

    end subroutine step_rk4

end module rk4

program nktg_demo
    use nktg_types
    use nktg_model
    use rk4
    implicit none

    ! ---------- Simulation setup ----------
    real(dp) :: t, tEnd, h
    integer :: nSteps, i

    ! Initial conditions
    real(dp) :: x0, v0, m0

    type(Params) :: par
    real(dp) :: y(3), ynew(3)
    real(dp) :: p, g1, g2

    ! Example parameters (adjust to your scenario)
    x0 = 0.0_dp ! m
    v0 = 5.0_dp ! m/s
    m0 = 10.0_dp ! kg

    par%k = 10.0_dp ! N/m, spring pulls toward 0
    par%c = 0.0_dp ! N*s/m, mild damping
    par%Fx = 0.0_dp ! no constant external force
    par%thrust = 0.0_dp ! constant thrust (N)
    par%mdot = 0.2_dp ! kg/s (fuel burn)

    h = 0.01_dp ! s, time step
    tEnd = 5.0_dp ! s, total time
    nSteps = int(tEnd/h)

    y = [x0, v0, m0]
    t = 0.0_dp

    ! ---------- CSV header ----------
    write(*,'(A)') 't,x,v,m,p,NKTg1,NKTg2'

    ! ---------- Time integration ----------
    do i = 0, nSteps
        p = momentum(y(2), y(3))
        g1 = nktg1(y(1), y(2), y(3)) / 10
        g2 = nktg2(par%mdot, y(2), y(3))
        write(*,'(F10.4,1A,ES16.8,1A,ES16.8,1A,ES16.8,1A,ES16.8,1A,ES16.8,1A,ES16.8)') &
        t, ",", y(1), ",", y(2), ",", y(3), ",", p, ",", g1, ",", g2

        if (i == nSteps) exit
        call step_rk4(t, y, h, par, rhs, ynew)
        y = ynew
        t = t + h

        ! Stop if fuel is exhausted (mass cannot drop below a threshold)
        if (y(3) <= 10.0_dp) then
            ! Recompute outputs for final truncated step print (optional)
            cycle  ! simulation continues but mass is clamped in rhs via max()
        end if
    end do
end program nktg_demo

Also one question, what’s the point of extracting xp and p*dm/dt as two separate quantities? what’s their significance?

Do let me know your thoughts.

2 Likes