Physics starter - Oscillator

MODEL "Physics starter - Oscillator"
DESCRIPTION "Underdamped motion with an analytical trajectory and first falling zero. Dimensionless working variables."
VERSION "1.0.0"
AUTHOR "Grid Team"
TAGS "physics", "starter", "oscillator"

# Inputs and governing equation x'' + 2 gamma x' + omega0^2 x = 0
input Damping = 0.1
  VALIDATE BETWEEN 0 AND 0.8 MESSAGE "The analytical reference is underdamped; use damping from 0 to 0.8."
input Frequency = 2
  VALIDATE BETWEEN 1 AND 4 MESSAGE "Use frequency from 1 to 4; the event interval brackets the first zero."
input EndTime = 10
  VALIDATE BETWEEN 0.1 AND 10 MESSAGE "Use an end time from 0.1 to 10."
NaturalFrequency = SQRT(Frequency^2-Damping^2)
output Trajectory = ODE_SOLVE((time,state) => [INDEX(state,1,2),-2*Damping*INDEX(state,1,2)-Frequency^2*INDEX(state,1,1)], [1,0], 0, EndTime, {samples:201,rtol:1e-9,atol:1e-12})
output FinalPosition = INDEX(Trajectory,201,2)
output ExactPosition = EXP(-Damping*EndTime)*(COS(NaturalFrequency*EndTime)+Damping/NaturalFrequency*SIN(NaturalFrequency*EndTime))
output AbsoluteError = ABS(FinalPosition-ExactPosition)
output FirstZero = ODE_EVENT_TIME((time,state) => [INDEX(state,1,2),-2*Damping*INDEX(state,1,2)-Frequency^2*INDEX(state,1,1)], [1,0], 0, 6, (time,state) => INDEX(state,1,1), {rtol:1e-9,atol:1e-12,event_direction:"falling"})
output ExactFirstZero = (PI()/2+ATAN(Damping/NaturalFrequency))/NaturalFrequency
A1 = "Position at final time"
B1 = FinalPosition
A2 = "Analytical position"
B2 = ExactPosition
A3 = "Absolute error"
B3 = AbsoluteError
A4 = "First falling zero"
B4 = FirstZero
A5 = "Analytical first zero"
B5 = ExactFirstZero

# Scientific checks and interpretation
output Comparison = MAKEARRAY(201, 3, (row, col) => col = 1 THEN INDEX(Trajectory, row, 1)
  ELSE col = 2 THEN INDEX(Trajectory, row, 2)
  ELSE EXP(-Damping * INDEX(Trajectory, row, 1)) * (
    COS(NaturalFrequency * INDEX(Trajectory, row, 1))
    + Damping / NaturalFrequency * SIN(NaturalFrequency * INDEX(Trajectory, row, 1))))
output MaximumTrajectoryError = MAX(MAP(SEQUENCE(201), (row) => ABS(INDEX(Comparison, row, 2) - INDEX(Comparison, row, 3))))

output Assumptions = "Underdamped motion with an analytical trajectory and first falling zero. Dimensionless working variables. Solver completion alone does not establish global accuracy or physical validity."
output ScenarioGuide = "Save a Baseline input scenario, change an admitted input, inspect numerical checks and the unsmoothed plot, then restore Baseline. Save and reopen a .gridoc archive. Plot selection also accepts the Comparison numeric array; choose its X and Y columns in the chart dialog."

Charts_1!type = "charts"
Charts_1!config = <toml>
  version = 1
  kind = "charts"
  title = "Oscillator"
  [layout]
  x_axis = "Time"
  y_axis = "Displacement"
  [[charts]]
  type = "line"
  smooth = false
  bind = "Comparison"
  x = "0"
  y = "1"
  label = "Numerical position"
  color = "#2563eb"
  [[charts]]
  type = "line"
  smooth = false
  bind = "Comparison"
  x = "0"
  y = "2"
  label = "Analytical position"
  color = "#d97706"
</toml>

END MODEL