MATLAB中deval函数的Python等效实现及疫情模型求解需求

Hey there! I see you're looking for a Python alternative to MATLAB's deval function to find when infected and recovered populations are equal in your SIR epidemic model. Let's break down two solid approaches you can use here:

Since Python doesn't have a direct equivalent to MATLAB's deval, we can use either interpolation + root-finding or direct root-solving with ODE evaluations to get the time point where i(t) = r(t). Here's how to implement both with your existing code:

Method 1: Interpolation + Root-Finding #

First, we'll solve the ODE once to get discrete state values, then create smooth continuous functions for i(t) and r(t) using interpolation. Finally, we'll find where these two functions intersect.

Code Implementation #

import numpy as np
from scipy.integrate import odeint
from scipy.interpolate import interp1d
from scipy.optimize import root_scalar
import matplotlib.pyplot as plt

# Define your SIR epidemic model
def epidemic_model(state, t):
    s, i, r = state
    d_s = -0.06 * s * i
    d_i = 0.06 * s * i - 7*i
    d_r = 7 * i
    return [d_s, d_i, d_r]

# Solve the ODE for discrete time points
t = np.arange(0, 1, 0.01)
init_state = [990, 10, 0]
state = odeint(epidemic_model, init_state, t)
s_vals, i_vals, r_vals = state.T

# Create smooth continuous interpolation functions for i(t) and r(t)
i_interp = interp1d(t, i_vals, kind='cubic')  # Cubic interpolation for smoothness
r_interp = interp1d(t, r_vals, kind='cubic')

# Define the function we want to find roots for: i(t) - r(t) = 0
def find_intersection(t_val):
    return i_interp(t_val) - r_interp(t_val)

# Use root-finding to locate the intersection point
result = root_scalar(find_intersection, bracket=[0, 1], method='bisect')

# Output and visualize the result
if result.converged:
    print(f"Time when infected = recovered: t = {result.root:.4f}")
    print(f"Population count at this time: i = r = {i_interp(result.root):.2f}")
    
    # Plot to verify
    t_fine = np.linspace(0, 1, 1000)
    plt.figure(figsize=(10,6))
    plt.plot(t_fine, i_interp(t_fine), label='Infected Population (i(t))')
    plt.plot(t_fine, r_interp(t_fine), label='Recovered Population (r(t))')
    plt.scatter(result.root, i_interp(result.root), color='red', zorder=5, label='i = r Intersection')
    plt.xlabel('Time')
    plt.ylabel('Population')
    plt.legend()
    plt.grid(True)
    plt.show()
else:
    print("Failed to find a valid intersection point.")

Method 2: Direct Root-Solving with ODE Evaluations #

This approach skips interpolation and directly solves for the time t where i(t) - r(t) = 0 by evaluating the ODE at candidate time points during the root-finding process. It's more accurate but slightly slower.

Code Implementation #

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

# Define your SIR epidemic model
def epidemic_model(state, t):
    s, i, r = state
    d_s = -0.06 * s * i
    d_i = 0.06 * s * i - 7*i
    d_r = 7 * i
    return [d_s, d_i, d_r]

init_state = [990, 10, 0]

# Define the target function: returns i(t) - r(t) for a given t
def target_func(t_val):
    t_arr = np.array([t_val])
    state = odeint(epidemic_model, init_state, t_arr)
    return state[0, 1] - state[0, 2]

# Initial guess for the root (we can estimate this from the discrete data)
initial_guess = 0.1
root = fsolve(target_func, initial_guess)[0]

# Output and visualize the result
print(f"Time when infected = recovered: t = {root:.4f}")
state_at_root = odeint(epidemic_model, init_state, [root])[0]
print(f"Population count at this time: i = r = {state_at_root[1]:.2f}")

# Plot to verify
t = np.arange(0, 1, 0.01)
state = odeint(epidemic_model, init_state, t)
s_vals, i_vals, r_vals = state.T

plt.figure(figsize=(10,6))
plt.plot(t, i_vals, label='Infected Population (i(t))')
plt.plot(t, r_vals, label='Recovered Population (r(t))')
plt.scatter(root, state_at_root[1], color='red', zorder=5, label='i = r Intersection')
plt.xlabel('Time')
plt.ylabel('Population')
plt.legend()
plt.grid(True)
plt.show()

Quick Comparison #

• Interpolation Method: Faster, since we only solve the ODE once. Great if you need to query multiple state values after finding the intersection.
• Direct Root-Solving: More precise, as it avoids interpolation errors. Better for scenarios where accuracy is critical.

内容的提问来源于stack exchange，提问作者Tommy Ofinger