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