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

1

! 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 :: 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 = 10.0_dp ! m
v0 = 0.0_dp ! m/s
m0 = 100.0_dp ! kg

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

h = 0.01_dp ! s, time step
tEnd = 50.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))
g2 = nktg2(par%mdot, y(2), y(3))
write(*,‘(F10.4,1H,ES16.8,1H,ES16.8,1H,ES16.8,1H,ES16.8,1H,ES16.8,1H,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

Quick explanation (to pin inside the post)

  • State: y = [x, v, m]

  • Sample dynamics (can be adapted depending on the problem):

    dx/dt = v
dv/dt = (Fx - k*x - c*v + thrust)/m
dm/dt = mdot    (mdot < 0 for fuel burn)

  • Quantities:

    p      = m * v
NKTg1  = x * p
NKTg2  = (dm/dt) * p

  • Integration: Runge–Kutta 4 (RK4); CSV output columns:

    t, x, v, m, p, NKTg1, NKTg2
2

Welcome to the forum, but that said, what is the purpose of your post? I see no question nor an announcement of sorts. So I am a trifle puzzled :slight_smile: .

1 Like
3

“Thanks for the warm welcome! The purpose of my post is to share a minimal working example of my project (the NKTg law with varying inertia, using RK4 and CSV output). I’d like to hear feedback from the community, whether in terms of code style, numerical methods, or possible improvements. I realize I should have made that clearer in the first post.”

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