// TypeScript usage example
import {printPollution} from 'diff-grok';
printPollution();
# Python reference calculation for comparisonimport numpy as np
import pandas as pd
from scipy.integrate import solve_ivp
# -------------------------
# 0. Solver method and tolerance
# -------------------------
METHOD = "Radau" # stiff solver
RTOL = 1e-6
ATOL = 1e-6
# -------------------------
# 1. Reaction constants
# -------------------------
k1 = 0.35
k2 = 26.6
k3 = 1.23e4
k4 = 8.6e-4
k5 = 8.2e-4
k6 = 1.5e4
k7 = 1.3e-4
k8 = 2.4e4
k9 = 1.65e4
k10 = 9e3
k11 = 0.022
k12 = 1.2e4
k13 = 1.88
k14 = 1.63e4
k15 = 4.8e6
k16 = 3.5e-4
k17 = 0.0175
k18 = 1e8
k19 = 4.44e11
k20 = 1240
k21 = 2.1
k22 = 5.78
k23 = 0.0474
k24 = 1780
k25 = 3.12
# -------------------------
# 2. Initial conditions
# -------------------------
y0 = [
0, # y1 NO2
0.2, # y2 NO
0, # y3 O3P
0.04, # y4 O3
0, # y5 HO2
0, # y6 OH
0.1, # y7 HCHO
0.3, # y8 CO
0.01, # y9 ALD
0, # y10 MEO2
0, # y11 C2O3
0, # y12 CO2
0, # y13 PAN
0, # y14 CH3O
0, # y15 HNO3
0, # y16 O1D
0.007, # y17 SO2
0, # y18 SO4
0, # y19 NO3
0 # y20 N2O5
]
# -------------------------
# 3. Define POLL ODE system
# -------------------------
def poll(t, y):
y = np.array(y)
r1 = k1*y[0]
r2 = k2*y[1]*y[3]
r3 = k3*y[4]*y[1]
r4 = k4*y[6]
r5 = k5*y[6]
r6 = k6*y[6]*y[5]
r7 = k7*y[8]
r8 = k8*y[8]*y[5]
r9 = k9*y[10]*y[1]
r10 = k10*y[10]*y[0]
r11 = k11*y[12]
r12 = k12*y[9]*y[1]
r13 = k13*y[13]
r14 = k14*y[0]*y[5]
r15 = k15*y[2]
r16 = k16*y[3]
r17 = k17*y[3]
r18 = k18*y[15]
r19 = k19*y[15]
r20 = k20*y[16]*y[5]
r21 = k21*y[18]
r22 = k22*y[18]
r23 = k23*y[0]*y[3]
r24 = k24*y[18]*y[0]
r25 = k25*y[19]
dy = np.zeros(20)
dy[0] = -(r1 + r10 + r14 + r23 + r24) + (r2 + r3 + r9 + r11 + r12 + r22 + r25)
dy[1] = -r2 - r3 - r9 - r12 + r1 + r21
dy[2] = -r15 + r1 + r17 + r19 + r22
dy[3] = -r2 - r16 - r17 - r23 + r15
dy[4] = -r3 + 2*r4 + r6 + r7 + r13 + r20
dy[5] = -r6 - r8 - r14 - r20 + r3 + 2*r18
dy[6] = -r4 - r5 - r6 + r13
dy[7] = r4 + r5 + r6 + r7
dy[8] = -r7 - r8
dy[9] = -r12 + r7 + r9
dy[10] = -r9 - r10 + r8 + r11
dy[11] = r9
dy[12] = -r11 + r10
dy[13] = -r13 + r12
dy[14] = r14
dy[15] = -r18 - r19 + r16
dy[16] = -r20
dy[17] = r20
dy[18] = -r21 - r22 - r24 + r23 + r25
dy[19] = -r25 + r24
return dy
# -------------------------
# 4. Time array
# -------------------------
t_start = 0
t_end = 60 # final time
t_step = 0.002
n_points = int(np.ceil((t_end - t_start)/t_step)) + 1
t_eval = np.linspace(t_start, t_end, n_points)
# -------------------------
# 5. Solve ODE
# -------------------------
sol = solve_ivp(
poll,
(t_start, t_end),
y0,
method=METHOD,
t_eval=t_eval,
rtol=RTOL,
atol=ATOL
)
print("Success:", sol.success)
print("Message:", sol.message)
# -------------------------
# 6. Save solution to CSV
# -------------------------
col_names = ['t'] + [f'y{i+1}' for i in range(20)]
data = np.column_stack((sol.t, sol.y.T))
df = pd.DataFrame(data, columns=col_names)
filename = f"POLL-using-{METHOD}.csv"
df.to_csv(filename, index=False)
print(f"Solution saved to '{filename}'")
Print solution of the POLL benchmark problem