Python provides powerful tools for solving ordinary differential equations. Let's start with a simple first-order ODE:
Example 1: Solve dy/dt = -y + t with initial condition y(0) = 0
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Define the differential equation dy/dt = -y + t
def dydt(y, t):
return -y + t
# Initial condition
y0 = 0
# Time points
t = np.linspace(0, 5, 100)
# Solve ODE
y = odeint(dydt, y0, t)
# Plot the solution
plt.figure(figsize=(8, 5))
plt.plot(t, y)
plt.xlabel('Time t')
plt.ylabel('y(t)')
plt.title('Solution of dy/dt = -y + t')
plt.grid(True)
plt.show()
# Print the analytical solution for comparison
print("Numerical solution at t=5:", y[-1][0])
print("Analytical solution at t=5:", 5 - 1 + np.exp(-5))Direction fields (slope fields) help visualize the behavior of differential equations without solving them explicitly.
import numpy as np
import matplotlib.pyplot as plt
# Create a grid of points
x = np.linspace(-3, 3, 20)
y = np.linspace(-3, 3, 20)
X, Y = np.meshgrid(x, y)
# Define the differential equation dy/dx = x - y
def dydx(x, y):
return x - y
# Calculate the slope at each point
U = np.ones_like(X)
V = dydx(X, Y)
# Normalize the arrows
N = np.sqrt(U**2 + V**2)
U, V = U/N, V/N
# Create the direction field plot
plt.figure(figsize=(8, 6))
plt.quiver(X, Y, U, V, alpha=0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Direction Field for dy/dx = x - y')
plt.grid(True)
plt.axis('equal')
plt.show()A separable equation can be written as dy/dx = f(x)g(y). Here's how to solve them numerically:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Example: dy/dx = y*cos(x), which separates to dy/y = cos(x)dx
def dydx(y, x):
return y * np.cos(x)
# Initial condition: y(0) = 1
y0 = 1
x = np.linspace(0, 2*np.pi, 100)
# Numerical solution
y_numerical = odeint(dydx, y0, x)
# Analytical solution: y = e^(sin(x))
y_analytical = np.exp(np.sin(x))
# Plot both solutions
plt.figure(figsize=(10, 6))
plt.plot(x, y_numerical, 'b-', label='Numerical Solution', linewidth=2)
plt.plot(x, y_analytical, 'r--', label='Analytical Solution', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Separable Equation: dy/dx = y*cos(x)')
plt.legend()
plt.grid(True)
plt.show()
print("Max error:", np.max(np.abs(y_numerical.flatten() - y_analytical)))Some equations that aren't initially separable can be transformed into separable form using substitutions.
An equation dy/dx = f(y/x) is homogeneous. Use the substitution v = y/x to make it separable.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Example: dy/dx = (x + y)/(x - y)
# This is homogeneous since it can be written as dy/dx = (1 + y/x)/(1 - y/x)
# Using substitution v = y/x, we get x*dv/dx = (1 + v)/(1 - v) - v
def dvdx(v, x):
return ((1 + v)/(1 - v) - v)/x
# Initial condition: y(1) = 0, so v(1) = 0
v0 = 0
x = np.linspace(1, 5, 100)
# Solve for v
v = odeint(dvdx, v0, x)
# Convert back to y
y = v * x.reshape(-1, 1)
# Also solve the original equation directly for comparison
def dydx_original(y, x):
return (x + y)/(x - y + 1e-10) # Add small value to avoid division by zero
y0_original = 0
y_direct = odeint(dydx_original, y0_original, x)
# Plot results
plt.figure(figsize=(10, 6))
plt.plot(x, y, 'b-', label='Using substitution v = y/x', linewidth=2)
plt.plot(x, y_direct, 'r--', label='Direct solution', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Homogeneous Equation: dy/dx = (x + y)/(x - y)')
plt.legend()
plt.grid(True)
plt.show()
# Verify the substitution worked
print(f"Max difference between methods: {np.max(np.abs(y.flatten() - y_direct.flatten())):.6f}")Example 2: Equations of the form dy/dx = f(ax + by + c)
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Example: dy/dx = (2x + y - 1)²
# Use substitution u = 2x + y - 1
# From u = 2x + y - 1, we get du/dx = 2 + dy/dx
# So du/dx = 2 + u², which gives du/(2 + u²) = dx
def dudx(u, x):
return 2 + u**2
# Initial condition: y(0) = 1, so u(0) = 2(0) + 1 - 1 = 0
u0 = 0
x = np.linspace(0, 1, 100)
# Solve for u
u = odeint(dudx, u0, x)
# Convert back to y: y = u - 2x + 1
y = u - 2*x.reshape(-1, 1) + 1
# Plot the solution
plt.figure(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.plot(x, y, 'b-', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of dy/dx = (2x + y - 1)²')
plt.grid(True)
# Plot the substitution variable u
plt.subplot(1, 2, 2)
plt.plot(x, u, 'r-', linewidth=2)
plt.xlabel('x')
plt.ylabel('u = 2x + y - 1')
plt.title('Substitution Variable')
plt.grid(True)
plt.tight_layout()
plt.show()An equation M(x,y)dx + N(x,y)dy = 0 is exact if ∂M/∂y = ∂N/∂x. When exact, there exists a function F(x,y) such that dF = M dx + N dy.
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
# Example: (2x + 3y)dx + (3x + 4y)dy = 0
# Check if exact: ∂M/∂y = 3, ∂N/∂x = 3 ✓
# The solution is F(x,y) = x² + 3xy + 2y² = C
# Create contour plot of the solution curves
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
X, Y = np.meshgrid(x, y)
# Function F(x,y)
F = X**2 + 3*X*Y + 2*Y**2
# Plot level curves
plt.figure(figsize=(10, 8))
contours = plt.contour(X, Y, F, levels=20, colors='blue')
plt.clabel(contours, inline=True, fontsize=8)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Exact Equation: (2x + 3y)dx + (3x + 4y)dy = 0')
plt.grid(True)
plt.axis('equal')
# Add direction field
x_field = np.linspace(-3, 3, 15)
y_field = np.linspace(-3, 3, 15)
X_field, Y_field = np.meshgrid(x_field, y_field)
# From M dx + N dy = 0, we get dy/dx = -M/N
M = 2*X_field + 3*Y_field
N = 3*X_field + 4*Y_field
dydx = -M/(N + 1e-10)
# Normalize for direction field
U = np.ones_like(X_field)
V = dydx
N_norm = np.sqrt(U**2 + V**2)
U, V = U/N_norm, V/N_norm
plt.quiver(X_field, Y_field, U, V, alpha=0.5, color='red')
plt.show()
# Verify exactness numerically
def verify_exact(x, y):
h = 1e-6
# M = 2x + 3y, N = 3x + 4y
dM_dy = ((2*x + 3*(y+h)) - (2*x + 3*y))/h
dN_dx = ((3*(x+h) + 4*y) - (3*x + 4*y))/h
return dM_dy, dN_dx
x_test, y_test = 1.0, 2.0
dM_dy, dN_dx = verify_exact(x_test, y_test)
print(f"At point ({x_test}, {y_test}):")
print(f"∂M/∂y = {dM_dy:.6f}")
print(f"∂N/∂x = {dN_dx:.6f}")
print(f"Difference: {abs(dM_dy - dN_dx):.10f}")Finding the Solution Function:
import sympy as sp
# Define symbolic variables
x, y = sp.symbols('x y')
# Define M and N
M = 2*x + 3*y
N = 3*x + 4*y
# Check exactness
dM_dy = sp.diff(M, y)
dN_dx = sp.diff(N, x)
print(f"∂M/∂y = {dM_dy}")
print(f"∂N/∂x = {dN_dx}")
print(f"Is exact? {dM_dy == dN_dx}")
# Find F(x,y) by integrating M with respect to x
F = sp.integrate(M, x)
print(f"\nF(x,y) = ∫M dx = {F} + g(y)")
# Find g(y) by differentiating F with respect to y and comparing with N
dF_dy = sp.diff(F, y)
print(f"∂F/∂y = {dF_dy} + g'(y)")
print(f"This should equal N = {N}")
# So g'(y) = N - ∂F/∂y
g_prime = N - dF_dy
print(f"\ng'(y) = {g_prime}")
# Integrate to find g(y)
g = sp.integrate(g_prime, y)
print(f"g(y) = {g}")
# Complete solution
F_complete = F + g
print(f"\nComplete solution: F(x,y) = {F_complete} = C")When an equation M dx + N dy = 0 is not exact, we can sometimes find an integrating factor μ(x,y) such that μM dx + μN dy = 0 is exact.
Common types of integrating factors:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# Example: (3x²y + 2xy)dx + (x³ + x²)dy = 0
# This is NOT exact, but has an integrating factor μ = 1/x²
# After multiplying by μ = 1/x², we get:
# (3y + 2y/x)dx + (x + 1)dy = 0
# Verify this is now exact
def verify_exactness():
# M_new = 3y + 2y/x, N_new = x + 1
# ∂M_new/∂y = 3 + 2/x
# ∂N_new/∂x = 1
# Still not exact with this example, let's use a better one
pass
# Better example: (y)dx + (2x - ye^y)dy = 0
# This has integrating factor μ = e^(-y)
# Original equation
def plot_original():
x = np.linspace(0.1, 5, 20)
y = np.linspace(0.1, 3, 20)
X, Y = np.meshgrid(x, y)
# From y dx + (2x - ye^y) dy = 0
# dy/dx = -y/(2x - ye^y)
dydx = -Y/(2*X - Y*np.exp(Y))
# Direction field
U = np.ones_like(X)
V = dydx
N = np.sqrt(U**2 + V**2)
U, V = U/N, V/N
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.quiver(X, Y, U, V, alpha=0.6)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Original (Non-exact) Equation')
plt.grid(True)
# After multiplying by μ = e^(-y)
# We get: ye^(-y)dx + (2xe^(-y) - y)dy = 0
# This is exact with solution F = xye^(-y) + ye^(-y) - e^(-y) = C
plt.subplot(1, 2, 2)
# Plot solution curves
F = X*Y*np.exp(-Y) + Y*np.exp(-Y) - np.exp(-Y)
contours = plt.contour(X, Y, F, levels=15)
plt.clabel(contours, inline=True, fontsize=8)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution Curves (After Integrating Factor)')
plt.grid(True)
plt.tight_layout()
plt.show()
plot_original()
# Demonstrate finding an integrating factor
print("Example: Finding an integrating factor")
print("Equation: y dx + (2x - ye^y) dy = 0")
print("M = y, N = 2x - ye^y")
print("\n∂M/∂y = 1")
print("∂N/∂x = 2")
print("Not exact since ∂M/∂y ≠ ∂N/∂x")
print("\nTrying μ = μ(y):")
print("(∂N/∂x - ∂M/∂y)/M = (2 - 1)/y = 1/y")
print("This is a function of y only!")
print("\nSolving dμ/dy = -μ/y gives μ = e^(-y)")
print("\nVerification: After multiplying by e^(-y), the equation becomes exact.")For the equation (2y² - 3x)dx + 2xy dy = 0, determine if it has an integrating factor of the form μ(x).
Linear first-order ODEs have the form dy/dx + P(x)y = Q(x). Bernoulli equations have the form dy/dx + P(x)y = Q(x)y^n.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Example: dy/dx + 2y/x = x²
# This is linear with P(x) = 2/x and Q(x) = x²
# The integrating factor is μ = e^(∫P(x)dx) = e^(2ln|x|) = x²
def solve_linear_ode(x_vals, y0, P_func, Q_func):
"""Solve linear ODE dy/dx + P(x)y = Q(x)"""
def dydx(y, x):
return Q_func(x) - P_func(x) * y
return odeint(dydx, y0, x_vals)
# Define P(x) and Q(x)
P = lambda x: 2/x
Q = lambda x: x**2
# Initial condition
y0 = 1
x = np.linspace(0.1, 3, 100)
# Numerical solution
y_numerical = solve_linear_ode(x, y0, P, Q)
# Analytical solution: y = (x³/5 + C/x²)
# With y(0.1) = 1, we can find C
C = (1 - 0.1**3/5) * 0.1**2
y_analytical = x**3/5 + C/x**2
# Plot comparison
plt.figure(figsize=(10, 6))
plt.plot(x, y_numerical, 'b-', label='Numerical Solution', linewidth=2)
plt.plot(x, y_analytical, 'r--', label='Analytical Solution', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Linear ODE: dy/dx + 2y/x = x²')
plt.legend()
plt.grid(True)
plt.show()
print(f"Max error: {np.max(np.abs(y_numerical.flatten() - y_analytical)):.6e}")import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Example: dy/dx - y/x = x²y²
# This is Bernoulli with P(x) = -1/x, Q(x) = x², n = 2
# Substitution: v = y^(1-n) = y^(-1)
# This transforms to: dv/dx + (1-n)P(x)v = (1-n)Q(x)
# For our example: dv/dx + v/x = -x²
def solve_bernoulli(x_vals, y0, P_func, Q_func, n):
"""Solve Bernoulli equation dy/dx + P(x)y = Q(x)y^n"""
if n == 1:
# This is just a linear equation
def dydx(y, x):
return Q_func(x)*y - P_func(x)*y
else:
# Transform to linear equation for v = y^(1-n)
v0 = y0**(1-n)
def dvdx(v, x):
return (1-n)*(Q_func(x) - P_func(x)*v)
v = odeint(dvdx, v0, x_vals)
return v**(1/(1-n))
return odeint(dydx, y0, x_vals)
# Original Bernoulli equation
def dydx_bernoulli(y, x):
return y/x + x**2 * y**2
# Solve using transformation
P = lambda x: -1/x
Q = lambda x: x**2
n = 2
y0 = 0.5
x = np.linspace(0.5, 3, 100)
# Numerical solution of original equation
y_direct = odeint(dydx_bernoulli, y0, x)
# Solution using Bernoulli transformation
y_transformed = solve_bernoulli(x, y0, P, Q, n)
# Plot both solutions
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(x, y_direct, 'b-', label='Direct Numerical Solution', linewidth=2)
plt.plot(x, y_transformed, 'r--', label='Bernoulli Transformation', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Bernoulli Equation: dy/dx - y/x = x²y²')
plt.legend()
plt.grid(True)
# Also plot 1/y to show the linear behavior after transformation
plt.subplot(1, 2, 2)
plt.plot(x, 1/y_direct, 'b-', linewidth=2)
plt.xlabel('x')
plt.ylabel('v = 1/y')
plt.title('After Transformation (Linear)')
plt.grid(True)
plt.tight_layout()
plt.show()What substitution transforms the Bernoulli equation dy/dx + y = xy³ into a linear equation?
A Riccati equation has the form dy/dx = q₀(x) + q₁(x)y + q₂(x)y². If we know one particular solution y₁, we can find the general solution.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Example: dy/dx = 1 + x² - 2xy + y²
# This can be rewritten as dy/dx = (1 + x²) - 2xy + y²
# One particular solution is y₁ = x
# If y₁ is a known solution, use substitution y = y₁ + 1/v
# This transforms the Riccati equation to a linear equation in v
def riccati_ode(y, x):
"""Riccati equation: dy/dx = 1 + x² - 2xy + y²"""
return 1 + x**2 - 2*x*y + y**2
# Verify that y₁ = x is a solution
x_test = np.linspace(0, 2, 100)
y1 = x_test
dy1dx = np.ones_like(x_test) # d(x)/dx = 1
riccati_value = riccati_ode(y1, x_test)
print("Verification that y₁ = x is a particular solution:")
print(f"Max |dy₁/dx - f(x,y₁)|: {np.max(np.abs(dy1dx - riccati_value)):.6e}")
# Now solve using the transformation y = x + 1/v
# This gives: dv/dx + v(2x - 2x) = -1
# Simplifying: dv/dx = -1, so v = -x + C
# General solution: y = x + 1/(-x + C)
def general_solution(x, C):
return x + 1/(-x + C)
# Plot several solution curves
plt.figure(figsize=(10, 8))
x = np.linspace(-2, 2, 100)
# Different values of C
C_values = [-3, -2, -1, 0, 1, 2, 3]
colors = plt.cm.viridis(np.linspace(0, 1, len(C_values)))
for C, color in zip(C_values, colors):
# Avoid singularity at x = C
x_safe = x[np.abs(x - C) > 0.1]
y = general_solution(x_safe, C)
plt.plot(x_safe, y, color=color, linewidth=2, label=f'C = {C}')
# Highlight the particular solution y = x
plt.plot(x, x, 'r--', linewidth=3, label='Particular solution y = x')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Riccati Equation: dy/dx = 1 + x² - 2xy + y²')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(True)
plt.xlim(-2, 2)
plt.ylim(-5, 5)
plt.tight_layout()
plt.show()
# Solve numerically for comparison
y0_values = [-2, -1, 0, 1, 2]
x_num = np.linspace(0, 1.5, 100)
plt.figure(figsize=(10, 6))
for y0 in y0_values:
y_num = odeint(riccati_ode, y0, x_num)
plt.plot(x_num, y_num, linewidth=2, label=f'y(0) = {y0}')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Numerical Solutions of the Riccati Equation')
plt.legend()
plt.grid(True)
plt.show()To solve a Riccati equation:
The existence and uniqueness theorem tells us when an initial value problem dy/dx = f(x,y), y(x₀) = y₀ has a unique solution.
Picard-Lindelöf Theorem: If f(x,y) and ∂f/∂y are continuous in a rectangle containing (x₀,y₀), then the IVP has a unique solution in some interval around x₀.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# Example 1: Unique solution
# dy/dx = y²/³, y(0) = 0
# Here f(x,y) = y²/³ is not Lipschitz continuous at y = 0
def example_non_unique():
# This IVP has multiple solutions:
# y = 0 (trivial solution)
# y = (x/3)³ for x ≥ 0
# y = -(|x|/3)³ for x < 0
# And infinitely many others!
x = np.linspace(-2, 2, 200)
plt.figure(figsize=(12, 8))
# Plot some of the non-unique solutions
plt.plot(x, np.zeros_like(x), 'b-', linewidth=2, label='y = 0')
# Solution that "turns on" at x = 0
y1 = np.where(x >= 0, (x/3)**3, 0)
plt.plot(x, y1, 'r-', linewidth=2, label='Turns on at x=0')
# Solution that "turns on" at x = 0.5
y2 = np.where(x >= 0.5, ((x-0.5)/3)**3, 0)
plt.plot(x, y2, 'g-', linewidth=2, label='Turns on at x=0.5')
# Solution that "turns on" at x = -0.5
y3 = np.where(x >= -0.5, ((x+0.5)/3)**3, 0)
plt.plot(x, y3, 'm-', linewidth=2, label='Turns on at x=-0.5')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Non-Unique Solutions: dy/dx = y²/³, y(0) = 0')
plt.legend()
plt.grid(True)
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)
example_non_unique()
plt.show()
# Example 2: Finite-time blow-up
# dy/dx = y², y(0) = 1
def example_blow_up():
# Analytical solution: y = 1/(1-x)
# Solution blows up at x = 1
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# Left plot: solution near blow-up
x1 = np.linspace(0, 0.99, 100)
y1 = 1/(1 - x1)
ax1.plot(x1, y1, 'b-', linewidth=2)
ax1.axvline(x=1, color='r', linestyle='--', label='Singularity at x=1')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_title('Solution of dy/dx = y²')
ax1.set_ylim(0, 50)
ax1.grid(True)
ax1.legend()
# Right plot: direction field
x = np.linspace(-0.5, 1.5, 20)
y = np.linspace(-2, 10, 20)
X, Y = np.meshgrid(x, y)
U = np.ones_like(X)
V = Y**2
N = np.sqrt(U**2 + V**2)
U, V = U/N, V/N
ax2.quiver(X, Y, U, V, alpha=0.6)
# Plot several solution curves
for y0 in [0.5, 1, 2, -1, -2]:
if y0 > 0:
x_sol = np.linspace(0, 1-1/y0-0.01, 100)
else:
x_sol = np.linspace(0, 2, 100)
y_sol = odeint(lambda y, x: y**2, y0, x_sol)
ax2.plot(x_sol, y_sol, linewidth=2)
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('Direction Field and Solution Curves')
ax2.set_xlim(-0.5, 1.5)
ax2.set_ylim(-2, 10)
ax2.grid(True)
plt.tight_layout()
example_blow_up()
plt.show()
# Interactive demonstration of Lipschitz condition
print("\nLipschitz Condition Check:")
print("For dy/dx = f(x,y) to have unique solution, f must satisfy:")
print("|f(x,y₁) - f(x,y₂)| ≤ L|y₁ - y₂| for some constant L")
print("\nExample: f(x,y) = 3y + sin(x)")
print("|∂f/∂y| = |3| = 3, so L = 3")
print("This guarantees unique solution!")
print("\nCounter-example: f(x,y) = y²/³")
print("∂f/∂y = (2/3)y^(-1/3) → ∞ as y → 0")
print("Not Lipschitz at y = 0, so uniqueness fails!")First-order ODEs model many real-world phenomena. Let's explore some practical applications.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Logistic growth model: dP/dt = rP(1 - P/K)
# r = growth rate, K = carrying capacity
def logistic_growth(P, t, r, K):
return r * P * (1 - P/K)
# Parameters
r = 0.5 # growth rate
K = 1000 # carrying capacity
P0_values = [50, 200, 500, 1500] # different initial populations
t = np.linspace(0, 20, 100)
plt.figure(figsize=(12, 5))
# Left plot: solution curves
plt.subplot(1, 2, 1)
for P0 in P0_values:
P = odeint(logistic_growth, P0, t, args=(r, K))
plt.plot(t, P, linewidth=2, label=f'P₀ = {P0}')
plt.axhline(y=K, color='r', linestyle='--', label=f'Carrying capacity K = {K}')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Logistic Growth Model')
plt.legend()
plt.grid(True)
# Right plot: phase line
plt.subplot(1, 2, 2)
P = np.linspace(0, 1500, 100)
dPdt = logistic_growth(P, 0, r, K)
plt.plot(P, dPdt, 'b-', linewidth=2)
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=K, color='r', linestyle='--', label=f'K = {K}')
plt.xlabel('Population P')
plt.ylabel('dP/dt')
plt.title('Phase Line Analysis')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# First-order reaction: A → B
# Rate equation: dA/dt = -kA
def reaction_rate(A, t, k):
return -k * A
# Parameters
k_values = [0.1, 0.5, 1.0, 2.0] # reaction rate constants
A0 = 100 # initial concentration
t = np.linspace(0, 10, 100)
plt.figure(figsize=(10, 6))
for k in k_values:
A = odeint(reaction_rate, A0, t, args=(k,))
plt.plot(t, A, linewidth=2, label=f'k = {k}')
# Plot half-life marker
t_half = np.log(2)/k
plt.plot(t_half, A0/2, 'o', markersize=8)
plt.text(t_half, A0/2 + 5, f't½={t_half:.2f}', ha='center')
plt.xlabel('Time')
plt.ylabel('Concentration of A')
plt.title('First-Order Chemical Reaction A → B')
plt.legend()
plt.grid(True)
plt.show()
# Temperature-dependent reaction (Arrhenius equation)
def arrhenius_reaction():
# k(T) = A * exp(-Ea/RT)
A = 1e10 # pre-exponential factor
Ea = 50000 # activation energy (J/mol)
R = 8.314 # gas constant
T = np.linspace(273, 373, 100) # temperature range (K)
k = A * np.exp(-Ea/(R*T))
plt.figure(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.plot(T - 273.15, k, 'b-', linewidth=2)
plt.xlabel('Temperature (°C)')
plt.ylabel('Rate constant k')
plt.title('Arrhenius Equation')
plt.grid(True)
plt.subplot(1, 2, 2)
plt.plot(1000/T, np.log(k), 'r-', linewidth=2)
plt.xlabel('1000/T (K⁻¹)')
plt.ylabel('ln(k)')
plt.title('Arrhenius Plot')
plt.grid(True)
plt.tight_layout()
plt.show()
arrhenius_reaction()import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# Newton's law of cooling: dT/dt = -k(T - T_env)
# T = temperature, T_env = environment temperature
def cooling_law(T, t, k, T_env):
return -k * (T - T_env)
# Coffee cooling example
T0 = 90 # initial temperature (°C)
T_env = 20 # room temperature (°C)
k_values = [0.05, 0.1, 0.2] # cooling constants
t = np.linspace(0, 60, 200) # time in minutes
plt.figure(figsize=(12, 5))
# Left plot: temperature over time
plt.subplot(1, 2, 1)
for k in k_values:
T = odeint(cooling_law, T0, t, args=(k, T_env))
plt.plot(t, T, linewidth=2, label=f'k = {k}')
plt.axhline(y=T_env, color='r', linestyle='--', label='Room temp')
plt.axhline(y=50, color='g', linestyle=':', label='Drinkable temp')
plt.xlabel('Time (minutes)')
plt.ylabel('Temperature (°C)')
plt.title("Newton's Law of Cooling - Coffee Example")
plt.legend()
plt.grid(True)
# Right plot: with variable environment temperature
plt.subplot(1, 2, 2)
# Environment temperature changes (e.g., taking coffee outside)
def T_env_variable(t):
if t < 20:
return 20 # indoors
elif t < 40:
return 5 # outside in winter
else:
return 20 # back indoors
def cooling_variable_env(T, t):
return -0.1 * (T - T_env_variable(t))
T_var = odeint(cooling_variable_env, T0, t)
# Plot environment temperature
t_env = t
T_env_plot = [T_env_variable(ti) for ti in t_env]
plt.plot(t, T_var, 'b-', linewidth=2, label='Coffee temperature')
plt.plot(t_env, T_env_plot, 'r--', linewidth=2, label='Environment temperature')
plt.xlabel('Time (minutes)')
plt.ylabel('Temperature (°C)')
plt.title('Cooling with Variable Environment')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()Modify the population growth model to include harvesting at a constant rate h. The equation becomes: dP/dt = rP(1-P/K) - h