Scribd
Upload
8 views2 pages

Import Numpy As NP

The document contains a Python script that models a set of chemical reactions using differential equations. It utilizes the SciPy library to solve these equations and plots the concentrations of products V and W over time. The script defines parameters, initial conditions, and a stopping condition for the simulation.

Uploaded by

35 Nguyễn Văn Thoại
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd
8 views2 pages

Import Numpy As NP

The document contains a Python script that models a set of chemical reactions using differential equations. It utilizes the SciPy library to solve these equations and plots the concentrations of products V and W over time. The script defines parameters, initial conditions, and a stopping condition for the simulation.

Uploaded by

35 Nguyễn Văn Thoại
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd
You are on page 1/ 2

import numpy as np

from scipy.integrate import solve_ivp

import matplotlib.pyplot as plt

# Parameters

k1 = 0.1

k2 = 0.01

k3 = 0.2

# Differential equations

def reactions(t, y):

Ca, Cb, Cv, Cw = y

rv = k1 * Ca * Cb

rw = k2 * Ca**2 / (1 + k3 * Cv**0.5)

dCa_dt = -rv - 2*rw

dCb_dt = -rv

dCv_dt = rv

dCw_dt = rw

return [dCa_dt, dCb_dt, dCv_dt, dCw_dt]

# Initial conditions

Ca0 = 0.25

Cb0 = 0.75

Cv0 = 0.0

Cw0 = 0.0

y0 = [Ca0, Cb0, Cv0, Cw0]

# Time span (up to 90% A reacted)

# 90% A reacted means Ca = 0.1*Ca0 = 0.025

# Solve until Ca reaches 0.025

def stop_condition(t, y):


return y[0] - 0.025

stop_condition.terminal = True

t_span = (0, 1000) # Initial time span, will be adjusted by the solver

t_eval = np.linspace(0, 1000, 100) # Evaluate at these points

# Solve the ODEs

solution = solve_ivp(reactions, t_span, y0, dense_output=True, events=stop_condition,


t_eval=t_eval)

# Extract results

t = solution.t

Ca = solution.y[0]

Cb = solution.y[1]

Cv = solution.y[2]

Cw = solution.y[3]

# Plotting

plt.figure(figsize=(10, 6))

plt.plot(t, Cv, label='V')

plt.plot(t, Cw, label='W')

plt.xlabel('Time (min)')

plt.ylabel('Concentration (M)')

plt.title('Concentration of V and W vs. Time')

plt.legend()

plt.grid(True)

plt.show()

You might also like