<jupyter_start><jupyter_text>Praca domowa - <NAME><jupyter_code>import sympy
sympy.init_printing()
t, lambda3a, lambdaa12, N4, N12, N16, dN4, dN12, dN16, dt = sympy.symbols('t, lambda_3a, lambda_a12, N4, N12, N16, dN4, dN12, dN16, dt', real=True)

eqs = [
    sympy.Eq(dN4/dt, -3*lambda3a * N4 **3 - lambdaa12 * N4 * N12),
    sympy.Eq(dN12/dt, lambda3a * N4 **3 - lambdaa12 * N4 * N12),
    sympy.Eq(dN16/dt, lambdaa12 * N4 * N12)
]
eqs
m, rho = sympy.symbols('m, rho', real=True)
X4, X12, X16, dX4, dX12, dX16 = sympy.symbols('X4, X12, X16, dX4, dX12, dX16', real=True)
Xeqs = [
    sympy.Eq(X4, m/rho*4*N4),
    sympy.Eq(X12, m/rho*12*N12),
    sympy.Eq(X16, m/rho*16*N16),
]
Xeqs
subs = {X4: dX4, X12: dX12, X16: dX16, N4: dN4, N12: dN12, N16: dN16}
dXeqs = [eq.subs(subs) for eq in Xeqs]
dXeqs
full_conservation = [sympy.Eq(X4 + X12 + X16, 1), sympy.Eq(dX4 + dX12 + dX16, 0)]
full_conservation
all_eqs = eqs + Xeqs + dXeqs + full_conservation
all_eqs
X_all_eqs = [eq.subs(sympy.solve(Xeqs, [N4, N12, N16])).subs(sympy.solve(dXeqs, [dN4, dN12, dN16])) for eq in eqs] + [full_conservation[1]]
X_all_eqs
solutions = sympy.solve(X_all_eqs, [dX4, dX12, dX16])
dX12dX4 = solutions[dX12]/solutions[dX4]
dX12dX4
q = sympy.symbols('q', real=True)
dX12dX4_final = dX12dX4.subs({lambdaa12*m: q * lambda3a * rho}).simplify()
dX12dX4_final
fX12 = sympy.Function('X12')(X4)
diffeq = sympy.Eq(fX12.diff(X4), dX12dX4_final.subs(X12, fX12))
diffeq
dX16dX4 = solutions[dX16]/solutions[dX4]
dX16dX4
dX16dX4_final = dX16dX4.subs({lambdaa12*m: q * lambda3a * rho}).simplify()
dX16dX4_final
derivatives_func = sympy.lambdify((X4, X12, X16, q), [dX12dX4_final, dX16dX4_final])
derivatives_func(1, 0, 0, 1)
def f(X, X4, q):
    return derivatives_func(X4, *X, q)
f([0, 0], 1, 1)
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
X4 = np.linspace(1, 0, 1000)
q_list = np.logspace(-3, np.log10(2), 500)
results = []
# fig, (ax1, ax2) = plt.subplots(2, sharex=True, figsize=(10, 8))
# ax1.set_xlim(0, 1)
# ax2.set_xlim(0, 1)
# ax1.set_ylim(0, 1)
# ax2.set_ylim(0, 1)
for q in q_list:
    X = odeint(f, [0, 0], X4, args=(q,))
    X12, X16 = X.T
#     ax1.plot(X4, X12, label=f"q: {q:.1f}")
#     ax2.plot(X4, X16, label=f"q: {q:.1f}")
#     ax2.set_xlabel("X4")
#     ax1.set_ylabel("X12")
#     ax2.set_ylabel("X16")
#     plt.plot(X4, X16)
#     plt.legend()
    results.append(X[-1])
results = np.array(results)
X12, X16 = results.T
plt.figure(figsize=(10, 10))
plt.plot(q_list, X12, label="X12")
plt.plot(q_list, X16, label="X16")
plt.xlabel("q")
plt.xscale("log")
plt.ylabel("X")
plt.legend(loc='best')
plt.xlim(q_list.min(), q_list.max());
plt.grid()
plt.savefig("Reacts.png")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19045871559633026, 'fraction_numerical': 0.10678899082568807, 'mean_word_length': 4.714884696016772, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '503079', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 1549, 'n_tokens_neox': 1415, 'n_words': 263, 'path': 'NuclearReactsOnly.ipynb', 'repo_name': 'StanczakDominik/FUW-AstroI-notatki', 'size': '2735'}