x(t) and T(t) for adiabatic piston on the gas oscillations in the elevator (https://phys.pro/problems/878/)

878.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc, rcParams
from scipy.integrate import solve_ivp

from matplotlib.ticker import FuncFormatter
from matplotlib.ticker import MultipleLocator

# Plot setup {{{
from matplotlib import rc, rcParams
rcParams['font.size'] = 45
rcParams['text.usetex'] = True
rcParams['font.sans-serif'] = ['Computer Modern']
rc('text.latex', preamble=r'\usepackage[T2A]{fontenc} \usepackage[utf8]{inputenc} \usepackage{bm}')

rc('axes', linewidth=1.8)
rc('xtick.major', size=4, width=2.0)
rc('ytick.major', size=4, width=2.0)
fig, ax = plt.subplots(figsize=(13,5))
fig.subplots_adjust(left=.20, bottom=.26, right=.90, top=.95)

# Set comma as a decimal separator (bad)
import locale
locale.setlocale(locale.LC_NUMERIC, 'ru_RU')
rcParams['axes.formatter.use_locale'] = True
def format_decimal(value, tick_number):
    value = round(value, 2) # Fixes 0,60000000000001 instead of 0,6
    return '$' + str(value).replace('.', '{,}') + '$'
ax.yaxis.set_major_formatter(FuncFormatter(format_decimal))
ax.xaxis.set_major_formatter(FuncFormatter(format_decimal))
#}}}

m = 1       # kg
g = 10      # m/s²
p0 = 1e5    # Pa
S = 0.01    # m²
a = 5       # m/s²
x1 = 1      # m
γ = 5/3     #
nu = 0.4    # mol
R = 8.31    # J / (mol·K)
F = 0.21    # N (Friction force)

p1 = m*g/S + p0
p2 = m*(g+a)/S + p0
x2 = x1*(p1/p2)**(1/γ)

print(x1, x2)

def f(t, y):
    x, v = y
    dxdt = v
    dvdt = (m*g + p0*S)/m * ((x1/x)**γ - 1) - a - F/m*np.sign(v)
    return [dxdt, dvdt]

t_span = (0, 1)
y0 = [x1, 0]
sol = solve_ivp(f, t_span, y0, t_eval=np.arange(0, 1, 1e-5), method='RK45')

time = sol.t
position = sol.y[0]
velocity = sol.y[1]
temperature = (m*g + p0*S) / (nu*R) * x1**γ * position**(1-γ)
ax.set_xlabel(r'$t, \textrm{с}$')
ax.grid(lw=1.8, color='k', ls='--')

## Uncomment to plot x(t)
#ax.plot(time, position*100, 'k', lw=3)
#ax.set_ylabel(r'$x, \textrm{см}$')
#plt.axhline(y=x1*100)
#ax.text(np.amax(time)*1.07, x1*100, '$x_1$', va='center', clip_on=False)
#plt.axhline(y=x2*100)
#ax.text(np.amax(time)*1.07, x2*100, '$x_2$', va='center', clip_on=False)

# Plot temperature T(t)
ax.plot(time, temperature, 'k', lw=3)
ax.set_ylabel(r'$T, \textrm{K}$')

plt.show()