Chapter 2: Maxwell Relations and Thermodynamic Identities

import sympy as sp import numpy as np import matplotlib.pyplot as plt # Define symbols T, S, P, V, mu, N = sp.symbols('T S P V mu N', real=True, positive=True) # Four thermodynamic potentials U = sp.Function('U')(S, V, N) F = sp.Function('F')(T, V, N) H = sp.Function('H')(S, P, N) G = sp.Function('G')(T, P, N) # Function to derive a Maxwell relation def derive_maxwell_relation(potential, var1, var2, potential_name): """Derive a Maxwell relation""" # First partial derivatives first_deriv_1 = sp.diff(potential, var1) first_deriv_2 = sp.diff(potential, var2) # Second partial derivatives (with the order swapped) second_deriv_12 = sp.diff(first_deriv_1, var2) second_deriv_21 = sp.diff(first_deriv_2, var1) # Maxwell relation maxwell_eq = sp.Eq(second_deriv_12, second_deriv_21) return maxwell_eq # Derive the four Maxwell relations (at constant particle number) print("=== Deriving the Maxwell Relations ===\n") # 1. From the internal energy U(S, V) print("1. Internal energy U(S, V):") print(" dU = T dS - P dV") print(" ∂T/∂V|_S = ∂(-P)/∂S|_V") print(" Maxwell relation: (∂T/∂V)_S = -(∂P/∂S)_V") print() # 2. From the Helmholtz free energy F(T, V) print("2. Helmholtz free energy F(T, V):") print(" dF = -S dT - P dV") print(" ∂(-S)/∂V|_T = ∂(-P)/∂T|_V") print(" Maxwell relation: (∂S/∂V)_T = (∂P/∂T)_V") print() # 3. From the enthalpy H(S, P) print("3. Enthalpy H(S, P):") print(" dH = T dS + V dP") print(" ∂T/∂P|_S = ∂V/∂S|_P") print(" Maxwell relation: (∂T/∂P)_S = (∂V/∂S)_P") print() # 4. From the Gibbs free energy G(T, P) print("4. Gibbs free energy G(T, P):") print(" dG = -S dT + V dP") print(" ∂(-S)/∂P|_T = ∂V/∂T|_P") print(" Maxwell relation: (∂S/∂P)_T = -(∂V/∂T)_P") print() # Summary table maxwell_relations = [ ("U(S,V)", "(∂T/∂V)_S", "-(∂P/∂S)_V"), ("F(T,V)", "(∂S/∂V)_T", "(∂P/∂T)_V"), ("H(S,P)", "(∂T/∂P)_S", "(∂V/∂S)_P"), ("G(T,P)", "(∂S/∂P)_T", "-(∂V/∂T)_P") ] print("=== List of Maxwell Relations ===") print(f"{'Potential':<12} {'LHS':<15} {'RHS':<15}") print("-" * 50) for pot, lhs, rhs in maxwell_relations: print(f"{pot:<12} {lhs:<15} = {rhs:<15}") # Practical example: verification with an ideal gas print("\n=== Verification with an Ideal Gas ===") print("Ideal gas: PV = NkT, U = (3/2)NkT (monatomic)") print("\nVerifying the Maxwell relation (∂S/∂V)_T = (∂P/∂T)_V:") print(" LHS: (∂S/∂V)_T = Nk/V") print(" RHS: (∂P/∂T)_V = Nk/V") print(" → They match!")

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import fsolve # van der Waals equation of state # (P + a/V²)(V - b) = RT def van_der_waals_pressure(V, T, a, b, R): """van der Waals pressure""" return R * T / (V - b) - a / V**2 def compute_heat_capacity_difference(V, T, a, b, R): """Compute C_P - C_V""" # (∂P/∂T)_V dP_dT_V = R / (V - b) # (∂P/∂V)_T dP_dV_T = -R * T / (V - b)**2 + 2 * a / V**3 # C_P - C_V = -T (∂P/∂T)_V² / (∂P/∂V)_T if dP_dV_T != 0: diff = -T * dP_dT_V**2 / dP_dV_T else: diff = np.nan return diff # van der Waals parameters for Ar (argon) R = 8.314 # J/(mol·K) a = 0.1355 # Pa·m⁶/mol² (1.355 bar·L²/mol²) b = 3.201e-5 # m³/mol (0.03201 L/mol) # Temperature range T_range = np.linspace(100, 500, 50) # K V_fixed = 1e-3 # 1 L/mol = 1e-3 m³/mol # Compute C_P - C_V Cp_minus_Cv_vdw = [] Cp_minus_Cv_ideal = [] for T in T_range: # van der Waals gas diff_vdw = compute_heat_capacity_difference(V_fixed, T, a, b, R) Cp_minus_Cv_vdw.append(diff_vdw) # Ideal gas (C_P - C_V = R) Cp_minus_Cv_ideal.append(R) # Visualization fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Temperature dependence of C_P - C_V ax1 = axes[0] ax1.plot(T_range, Cp_minus_Cv_vdw, 'b-', linewidth=2, label='van der Waals') ax1.plot(T_range, Cp_minus_Cv_ideal, 'r--', linewidth=2, label='Ideal gas (= R)') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('C_P - C_V (J/(mol·K))') ax1.set_title('Temperature Dependence of C_P - C_V (Ar, V = 1 L/mol)') ax1.legend() ax1.grid(True, alpha=0.3) # Volume dependence V_range = np.linspace(5e-5, 5e-3, 100) # m³/mol T_fixed = 300 # K Cp_minus_Cv_vs_V = [] for V in V_range: diff = compute_heat_capacity_difference(V, T_fixed, a, b, R) Cp_minus_Cv_vs_V.append(diff) ax2 = axes[1] ax2.plot(V_range * 1000, Cp_minus_Cv_vs_V, 'g-', linewidth=2) ax2.axhline(R, color='r', linestyle='--', linewidth=2, label='Ideal gas') ax2.set_xlabel('Molar volume (L/mol)') ax2.set_ylabel('C_P - C_V (J/(mol·K))') ax2.set_title(f'Volume Dependence of C_P - C_V (Ar, T = {T_fixed} K)') ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('thermo_heat_capacity_difference.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("=== Heat Capacity Relation (Ar at 300 K, 1 L/mol) ===") T = 300 V = 1e-3 diff_vdw = compute_heat_capacity_difference(V, T, a, b, R) print(f"van der Waals: C_P - C_V = {diff_vdw:.4f} J/(mol·K)") print(f"Ideal gas: C_P - C_V = {R:.4f} J/(mol·K)") print(f"Relative error: {abs(diff_vdw - R) / R * 100:.2f}%") print("\nIn the low-density limit (V → ∞) the gas approaches ideal behavior")

import numpy as np import matplotlib.pyplot as plt # Compute the isothermal compressibility def isothermal_compressibility(V, T, a, b, R): """Isothermal compressibility κ_T = -1/V (∂V/∂P)_T""" # Compute (∂P/∂V)_T and take the reciprocal dP_dV_T = -R * T / (V - b)**2 + 2 * a / V**3 if dP_dV_T != 0: kappa_T = -1 / (V * dP_dV_T) else: kappa_T = np.inf return kappa_T def volumetric_expansion(V, T, a, b, R): """Volumetric expansion coefficient α = 1/V (∂V/∂T)_P""" # Maxwell relation: (∂V/∂T)_P = -(∂S/∂P)_T # For simplicity, however, compute from (∂P/∂T)_V and (∂P/∂V)_T dP_dT_V = R / (V - b) dP_dV_T = -R * T / (V - b)**2 + 2 * a / V**3 if dP_dV_T != 0: dV_dT_P = -dP_dT_V / dP_dV_T alpha = dV_dT_P / V else: alpha = np.inf return alpha # Parameters for Ar R = 8.314 a = 0.1355 b = 3.201e-5 # Temperature and volume ranges T_range = np.linspace(150, 500, 100) V_range = np.linspace(1e-4, 5e-3, 100) # Temperature dependence (fixed V) V_fixed = 1e-3 kappa_T_vs_T = [] alpha_vs_T = [] for T in T_range: kappa_T = isothermal_compressibility(V_fixed, T, a, b, R) alpha = volumetric_expansion(V_fixed, T, a, b, R) kappa_T_vs_T.append(kappa_T) alpha_vs_T.append(alpha) # Volume dependence (fixed T) T_fixed = 300 kappa_T_vs_V = [] alpha_vs_V = [] for V in V_range: kappa_T = isothermal_compressibility(V, T_fixed, a, b, R) alpha = volumetric_expansion(V, T_fixed, a, b, R) kappa_T_vs_V.append(kappa_T) alpha_vs_V.append(alpha) # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Temperature dependence of the isothermal compressibility ax1 = axes[0, 0] ax1.plot(T_range, np.array(kappa_T_vs_T) * 1e9, 'b-', linewidth=2) ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('κ_T (GPa⁻¹)') ax1.set_title(f'Temperature Dependence of κ_T (V = {V_fixed*1000:.1f} L/mol)') ax1.grid(True, alpha=0.3) # Temperature dependence of the volumetric expansion coefficient ax2 = axes[0, 1] ax2.plot(T_range, np.array(alpha_vs_T) * 1e3, 'r-', linewidth=2) ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('α (10⁻³ K⁻¹)') ax2.set_title(f'Temperature Dependence of α (V = {V_fixed*1000:.1f} L/mol)') ax2.grid(True, alpha=0.3) # Volume dependence of the isothermal compressibility ax3 = axes[1, 0] ax3.plot(V_range * 1000, np.array(kappa_T_vs_V) * 1e9, 'g-', linewidth=2) ax3.set_xlabel('Molar volume (L/mol)') ax3.set_ylabel('κ_T (GPa⁻¹)') ax3.set_title(f'Volume Dependence of κ_T (T = {T_fixed} K)') ax3.grid(True, alpha=0.3) # Volume dependence of the volumetric expansion coefficient ax4 = axes[1, 1] ax4.plot(V_range * 1000, np.array(alpha_vs_V) * 1e3, 'm-', linewidth=2) ax4.set_xlabel('Molar volume (L/mol)') ax4.set_ylabel('α (10⁻³ K⁻¹)') ax4.set_title(f'Volume Dependence of α (T = {T_fixed} K)') ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('thermo_compressibility_expansion.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("=== Compressibility and Expansion Coefficient (Ar at 300 K, 1 L/mol) ===") T = 300 V = 1e-3 kappa_T = isothermal_compressibility(V, T, a, b, R) alpha = volumetric_expansion(V, T, a, b, R) print(f"Isothermal compressibility κ_T = {kappa_T*1e9:.4f} GPa⁻¹") print(f"Volumetric expansion coefficient α = {alpha*1e3:.4f} × 10⁻³ K⁻¹") # Comparison with an ideal gas kappa_T_ideal = V / (R * T) # Ideal gas: κ_T = 1/P alpha_ideal = 1 / T # Ideal gas: α = 1/T print(f"\nIdeal gas:") print(f" κ_T = {kappa_T_ideal*1e9:.4f} GPa⁻¹") print(f" α = {alpha_ideal*1e3:.4f} × 10⁻³ K⁻¹")

import numpy as np import sympy as sp # Define symbols T, P, V, S = sp.symbols('T P V S', real=True) def jacobian_2x2(u, v, x, y): """Compute a 2×2 Jacobian determinant""" J = sp.Matrix([ [sp.diff(u, x), sp.diff(u, y)], [sp.diff(v, x), sp.diff(v, y)] ]) return J.det() # Example: express (∂P/∂T)_V in terms of (∂S/∂V)_T # Derive the Maxwell relation: (∂S/∂V)_T = (∂P/∂T)_V print("=== Deriving Maxwell Relations with the Jacobian Method ===\n") # Treat S and V as functions of the variables S_func = sp.Function('S')(T, P) V_func = sp.Function('V')(T, P) # Properties of the Jacobian determinant print("1. Chain rule for Jacobians:") print(" ∂(S,V)/∂(T,P) · ∂(T,P)/∂(S,V) = 1") print() # Concrete example: from the Gibbs free energy print("2. Gibbs free energy G(T,P):") print(" dG = -S dT + V dP") print(" → S = -(∂G/∂T)_P, V = (∂G/∂P)_T") print() print("3. Commutativity of second partial derivatives:") print(" ∂²G/∂T∂P = ∂²G/∂P∂T") print(" → ∂S/∂P|_T = -∂V/∂T|_P") print(" This is a Maxwell relation") print() # Deriving practical identities print("=== Deriving Practical Identities ===\n") print("4. Express (∂U/∂V)_T in terms of measurable quantities:") print(" From dU = TdS - PdV,") print(" (∂U/∂V)_T = T(∂S/∂V)_T - P") print(" Using the Maxwell relation (∂S/∂V)_T = (∂P/∂T)_V,") print(" (∂U/∂V)_T = T(∂P/∂T)_V - P") print() # Numerical example: verify with a van der Waals gas from scipy.misc import derivative def U_vdw(V, T, a, b, R, n=1): """Internal energy of a van der Waals gas""" # U = (3/2)nRT - n²a/V (monatomic) return 1.5 * n * R * T - n**2 * a / V def P_vdw(V, T, a, b, R, n=1): """Pressure of a van der Waals gas""" return n * R * T / (V - n*b) - n**2 * a / V**2 # Parameters for Ar R = 8.314 a = 0.1355 b = 3.201e-5 n = 1 # 1 mol T = 300 V = 1e-3 # LHS: numerical derivative of (∂U/∂V)_T dU_dV_T = derivative(lambda v: U_vdw(v, T, a, b, R, n), V, dx=1e-8) # RHS: T(∂P/∂T)_V - P dP_dT_V = derivative(lambda t: P_vdw(V, t, a, b, R, n), T, dx=1e-6) P = P_vdw(V, T, a, b, R, n) right_hand_side = T * dP_dT_V - P print("5. Verification with a van der Waals gas (Ar, 300 K, 1 L/mol):") print(f" LHS (∂U/∂V)_T = {dU_dV_T:.4f} J/m³") print(f" RHS T(∂P/∂T)_V - P = {right_hand_side:.4f} J/m³") print(f" Relative error: {abs(dU_dV_T - right_hand_side)/abs(dU_dV_T)*100:.6f}%") print("\n → The identity is confirmed numerically!")

import numpy as np import matplotlib.pyplot as plt # Theoretical relation: κ_T / κ_S = C_P / C_V = γ def compute_gamma_ratio(V, T, a, b, R): """Compute γ = κ_T / κ_S = C_P / C_V""" # Isothermal compressibility dP_dV_T = -R * T / (V - b)**2 + 2 * a / V**3 kappa_T = -1 / (V * dP_dV_T) if dP_dV_T != 0 else np.inf # C_P - C_V dP_dT_V = R / (V - b) Cp_minus_Cv = -T * dP_dT_V**2 / dP_dV_T if dP_dV_T != 0 else np.inf # Correct relative to the monatomic ideal gas value C_V = (3/2)R Cv = 1.5 * R # Monatomic gas Cp = Cv + Cp_minus_Cv gamma = Cp / Cv # Adiabatic compressibility (from the theoretical relation) kappa_S = kappa_T / gamma return gamma, kappa_T, kappa_S, Cp, Cv # Parameters for Ar R = 8.314 a = 0.1355 b = 3.201e-5 # Compute γ over the temperature range T_range = np.linspace(150, 500, 100) V_fixed = 1e-3 gamma_values = [] kappa_T_values = [] kappa_S_values = [] Cp_values = [] Cv_values = [] for T in T_range: gamma, kappa_T, kappa_S, Cp, Cv = compute_gamma_ratio(V_fixed, T, a, b, R) gamma_values.append(gamma) kappa_T_values.append(kappa_T) kappa_S_values.append(kappa_S) Cp_values.append(Cp) Cv_values.append(Cv) # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # γ = C_P / C_V ax1 = axes[0, 0] ax1.plot(T_range, gamma_values, 'b-', linewidth=2) ax1.axhline(5/3, color='r', linestyle='--', linewidth=1.5, label='Ideal gas (5/3)') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('γ = C_P / C_V') ax1.set_title('Temperature Dependence of the Heat Capacity Ratio') ax1.legend() ax1.grid(True, alpha=0.3) # C_P and C_V ax2 = axes[0, 1] ax2.plot(T_range, Cp_values, 'r-', linewidth=2, label='C_P') ax2.plot(T_range, Cv_values, 'b-', linewidth=2, label='C_V') ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('Heat capacity (J/(mol·K))') ax2.set_title('Temperature Dependence of the Heat Capacities') ax2.legend() ax2.grid(True, alpha=0.3) # κ_T and κ_S ax3 = axes[1, 0] ax3.plot(T_range, np.array(kappa_T_values) * 1e9, 'g-', linewidth=2, label='κ_T') ax3.plot(T_range, np.array(kappa_S_values) * 1e9, 'm-', linewidth=2, label='κ_S') ax3.set_xlabel('Temperature (K)') ax3.set_ylabel('Compressibility (GPa⁻¹)') ax3.set_title('Temperature Dependence of the Compressibilities') ax3.legend() ax3.grid(True, alpha=0.3) # Verify the relation: κ_T / κ_S vs C_P / C_V ax4 = axes[1, 1] kappa_ratio = np.array(kappa_T_values) / np.array(kappa_S_values) Cp_Cv_ratio = np.array(Cp_values) / np.array(Cv_values) ax4.plot(T_range, kappa_ratio, 'b-', linewidth=2, label='κ_T / κ_S') ax4.plot(T_range, Cp_Cv_ratio, 'r--', linewidth=2, label='C_P / C_V') ax4.set_xlabel('Temperature (K)') ax4.set_ylabel('Ratio') ax4.set_title('Verification of κ_T / κ_S = C_P / C_V') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('thermo_kappa_gamma_relation.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("=== Verification of κ_T / κ_S = C_P / C_V (Ar at 300 K) ===") T = 300 V = 1e-3 gamma, kappa_T, kappa_S, Cp, Cv = compute_gamma_ratio(V, T, a, b, R) print(f"C_P = {Cp:.4f} J/(mol·K)") print(f"C_V = {Cv:.4f} J/(mol·K)") print(f"γ = C_P / C_V = {gamma:.6f}") print() print(f"κ_T = {kappa_T*1e9:.4f} GPa⁻¹") print(f"κ_S = {kappa_S*1e9:.4f} GPa⁻¹") print(f"κ_T / κ_S = {kappa_T/kappa_S:.6f}") print() print(f"Relative error: {abs(gamma - kappa_T/kappa_S)/gamma*100:.6f}%") print("\nFor a monatomic ideal gas, γ = 5/3 = 1.667") print(f"van der Waals gas: γ = {gamma:.4f}")

import numpy as np import matplotlib.pyplot as plt def check_stability(V, T, a, b, R): """Check the thermodynamic stability conditions""" # Check (∂P/∂V)_T < 0 dP_dV_T = -R * T / (V - b)**2 + 2 * a / V**3 # Check κ_T > 0 kappa_T = -1 / (V * dP_dV_T) if dP_dV_T != 0 else -np.inf # Check C_P > 0 (simplified: assume C_V > 0) Cv = 1.5 * R dP_dT_V = R / (V - b) Cp_minus_Cv = -T * dP_dT_V**2 / dP_dV_T if dP_dV_T != 0 else np.inf Cp = Cv + Cp_minus_Cv # Stability assessment stable = (dP_dV_T < 0) and (kappa_T > 0) and (Cp > 0) return { 'dP_dV_T': dP_dV_T, 'kappa_T': kappa_T, 'Cp': Cp, 'stable': stable } # Parameters for Ar R = 8.314 a = 0.1355 b = 3.201e-5 # Stability map in the T-V plane T_range = np.linspace(100, 500, 100) V_range = np.linspace(5e-5, 5e-3, 100) T_grid, V_grid = np.meshgrid(T_range, V_range) stability_map = np.zeros_like(T_grid) dP_dV_T_map = np.zeros_like(T_grid) for i in range(len(V_range)): for j in range(len(T_range)): V = V_range[i] T = T_range[j] result = check_stability(V, T, a, b, R) stability_map[i, j] = 1 if result['stable'] else 0 dP_dV_T_map[i, j] = result['dP_dV_T'] # Visualization fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Stability map ax1 = axes[0] c1 = ax1.contourf(T_grid, V_grid * 1000, stability_map, levels=[0, 0.5, 1], colors=['red', 'green'], alpha=0.5) ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('Molar volume (L/mol)') ax1.set_title('Thermodynamic Stability Map (green = stable, red = unstable)') # Map of (∂P/∂V)_T ax2 = axes[1] levels = np.linspace(-1e8, 1e8, 20) c2 = ax2.contourf(T_grid, V_grid * 1000, dP_dV_T_map, levels=levels, cmap='RdBu_r') ax2.contour(T_grid, V_grid * 1000, dP_dV_T_map, levels=[0], colors='black', linewidths=2) ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('Molar volume (L/mol)') ax2.set_title('(∂P/∂V)_T Map (negative = stable, positive = unstable)') plt.colorbar(c2, ax=ax2, label='(∂P/∂V)_T (Pa/m³)') plt.tight_layout() plt.savefig('thermo_stability_conditions.png', dpi=300, bbox_inches='tight') plt.show() # Compute the spinodal line ((∂P/∂V)_T = 0) print("=== Thermodynamic Stability Conditions ===\n") print("Spinodal line: the condition where (∂P/∂V)_T = 0") print("van der Waals: -RT/(V-b)² + 2a/V³ = 0") print("→ V_spinodal³ = 2a(V-b)² / (RT)") print() # Check at the critical point T_c = 8 * a / (27 * R * b) V_c = 3 * b print(f"Critical point (Ar):") print(f" T_c = {T_c:.2f} K") print(f" V_c = {V_c*1000:.4f} L/mol") result_c = check_stability(V_c, T_c, a, b, R) print(f" (∂P/∂V)_T = {result_c['dP_dV_T']:.2e} Pa/m³") print(f" → At the critical point, (∂P/∂V)_T = 0 (the stability boundary)")

import numpy as np import matplotlib.pyplot as plt from scipy.integrate import simpson # (Mock) experimental data: thermal expansion coefficient of Ar gas # α = 1/V (∂V/∂T)_P def generate_mock_data(): """Generate mock experimental data""" T_data = np.linspace(100, 400, 31) # 10 K intervals # Generated from the van der Waals model R = 8.314 a = 0.1355 b = 3.201e-5 P = 1e5 # 1 bar = 1e5 Pa alpha_data = [] for T in T_data: # Compute the volume at constant P (simplified: ideal-gas approximation) V = R * T / P # Thermal expansion coefficient dP_dT_V = R / (V - b) dP_dV_T = -R * T / (V - b)**2 + 2 * a / V**3 dV_dT_P = -dP_dT_V / dP_dV_T alpha = dV_dT_P / V alpha_data.append(alpha) return T_data, np.array(alpha_data) # Obtain the experimental data T_data, alpha_data = generate_mock_data() # Compute the entropy using the Maxwell relation: (∂S/∂P)_T = -(∂V/∂T)_P # ΔS = ∫(∂S/∂P)_T dP = -∫(∂V/∂T)_P dP = -∫V·α dP def compute_entropy_change(T, alpha, V, P_initial, P_final, n_points=100): """Entropy change accompanying a pressure change""" P_range = np.linspace(P_initial, P_final, n_points) # Volume at each pressure (ideal-gas approximation) R = 8.314 V_range = R * T / P_range # (∂S/∂P)_T = -V·α integrand = -V_range * alpha # Simpson integration delta_S = simpson(integrand, x=P_range) return delta_S # Compute the entropy change at each temperature P_initial = 1e5 # 1 bar P_final = 10e5 # 10 bar R = 8.314 entropy_changes = [] for i, T in enumerate(T_data): V = R * T / P_initial alpha = alpha_data[i] delta_S = compute_entropy_change(T, alpha, V, P_initial, P_final) entropy_changes.append(delta_S) # Visualization fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Thermal expansion coefficient data ax1 = axes[0] ax1.plot(T_data, alpha_data * 1e3, 'bo-', markersize=4, linewidth=2, label='Experimental data') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('α (10⁻³ K⁻¹)') ax1.set_title('Experimental Thermal Expansion Data') ax1.legend() ax1.grid(True, alpha=0.3) # Entropy change ax2 = axes[1] ax2.plot(T_data, entropy_changes, 'ro-', markersize=4, linewidth=2) ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('ΔS (J/(mol·K))') ax2.set_title(f'Entropy Change with Pressure\n({P_initial/1e5:.0f} → {P_final/1e5:.0f} bar)') ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('thermo_entropy_from_expansion.png', dpi=300, bbox_inches='tight') plt.show() # Display results print("=== Entropy Calculation from Experimental Data ===\n") print("Maxwell relation: (∂S/∂P)_T = -(∂V/∂T)_P = -V·α") print(f"Pressure change: {P_initial/1e5:.0f} bar → {P_final/1e5:.0f} bar\n") print(f"{'T (K)':<12} {'α (10⁻³ K⁻¹)':<18} {'ΔS (J/(mol·K))':<18}") print("-" * 50) for i in [0, 10, 20, 30]: print(f"{T_data[i]:<12.1f} {alpha_data[i]*1e3:<18.4f} {entropy_changes[i]:<18.4f}") print("\nComparison with theory (ideal gas):") T_ref = 300 delta_S_ideal = -R * np.log(P_final / P_initial) print(f" Ideal gas: ΔS = -R ln(P_f/P_i) = {delta_S_ideal:.4f} J/(mol·K)") print(f" Computed value (at {T_ref} K): {entropy_changes[20]:.4f} J/(mol·K)")