如何在scipy solve_ivp/odeint中输出并绘制时变/解依赖参数

Got it, let's tackle this problem step by step. The key issue here is that when you print b1 directly inside the vector field function, you're seeing output from every internal step that solve_ivp (or odeint) takes, not just the time points you care about. Here's how to fix this and get clean, time-aligned values for b1 and b2 to print and plot:

Core Approach #

Instead of printing parameters during the ODE solving process (which uses adaptive steps that don't match your output time points), we'll:

1. First solve the ODE to get state variables (x1, x2, etc.) at your desired time stamps
2. Recompute b1 and b2 using those solved state values—this ensures perfect alignment with your time points
3. Print and plot the aligned parameters alongside your solution

Modified Code for solve_ivp #

First, we'll adjust the vector field function to match solve_ivp's expected signature (t first, then state vector), then handle parameter calculation post-solve:

def vectorfield(t, w, p):
    """
    Defines the differential equations for the coupled spring-mass system.
    Args:
        t: Current time
        w: State vector [x1, dx1/dt, x2, dx2/dt]
        p: Fixed parameters [m1, m2, k1, k2, L1, L2]
    """
    x1, y1, x2, y2 = w
    m1, m2, k1, k2, L1, L2 = p
    
    # Calculate state-dependent friction coefficients
    b1 = 0.8 * x1
    b2 = 20 * x2 if x2 >= 1.0 else 1.0
    
    # Build the derivative vector
    derivatives = [
        y1,
        (-b1 * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
        y2,
        (-b2 * y2 - k2 * (x2 - x1 - L2)) / m2
    ]
    return derivatives

# Import required libraries
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# ----------------------
# Parameter & Setup
# ----------------------
# Fixed system parameters
m1 = 1.0
m2 = 1.5
k1 = 8.0
k2 = 40.0
L1 = 0.5
L2 = 1.0

# Initial conditions: [x1, dx1/dt, x2, dx2/dt]
w0 = [0.5, 0.0, 2.25, 0.0]

# Time settings: solve from 0 to stoptime, output numpoints evenly spaced points
stoptime = 10.0
numpoints = 101
t_eval = np.linspace(0, stoptime, numpoints)

# Pack fixed parameters
params = [m1, m2, k1, k2, L1, L2]

# ----------------------
# Solve the ODE
# ----------------------
sol = solve_ivp(
    vectorfield,
    t_span=[0, stoptime],
    y0=w0,
    args=(params,),
    t_eval=t_eval,  # Force output at our desired time points
    atol=1e-8,
    rtol=1e-6
)

# Extract solved states and time
t = sol.t
x1 = sol.y[0]
x2 = sol.y[2]

# ----------------------
# Calculate aligned b1 & b2
# ----------------------
# Compute b1 for each time point using solved x1
b1_vals = 0.8 * x1
# Compute b2 using vectorized condition (faster than loop)
b2_vals = np.where(x2 >= 1.0, 20 * x2, 1.0)

# ----------------------
# Print Time-Aligned Parameters
# ----------------------
print("Time (t) | b1 Value | b2 Value")
print("-------------------------------")
# Print first 10 points as an example (remove slice to print all)
for ti, b1i, b2i in zip(t[:10], b1_vals[:10], b2_vals[:10]):
    print(f"{ti:.2f}    | {b1i:.4f}   | {b2i:.4f}")

# ----------------------
# Plot Results
# ----------------------
plt.figure(figsize=(12, 8))

# Subplot 1: State variables (displacements and velocities)
plt.subplot(2, 1, 1)
plt.plot(t, x1, 'b', label=r'$x_1$ (Displacement)')
plt.plot(t, sol.y[1], 'g', label=r'$\dot{x_1}$ (Velocity)')
plt.plot(t, x2, 'r', label=r'$x_2$ (Displacement)')
plt.plot(t, sol.y[3], 'c', label=r'$\dot{x_2}$ (Velocity)')
plt.xlabel('Time (t)')
plt.grid(True)
plt.legend(fontsize=12)
plt.title('Coupled Spring-Mass System: State Variables')

# Subplot 2: Friction coefficients
plt.subplot(2, 1, 2)
plt.plot(t, b1_vals, 'k', label=r'$b_1 = 0.8x_1$')
plt.plot(t, b2_vals, 'm', label=r'$b_2 = \begin{cases}20x_2 & x_2 \geq 1.0 \\ 1.0 & x_2 < 1.0\end{cases}$')
plt.xlabel('Time (t)')
plt.grid(True)
plt.legend(fontsize=12)
plt.title('Time/State-Dependent Friction Coefficients')

plt.tight_layout()
plt.show()

Why This Works #

• solve_ivp uses adaptive step sizes to maintain accuracy, which means it calculates many more internal points than the ones you specify in t_eval. Printing inside the vector field function would show all these internal steps, leading to messy, unaligned output.
• By recomputing b1 and b2 after solving, we use exactly the state values at your desired time points—so every parameter value maps perfectly to a time stamp you care about.

If you ever need to track parameters during the internal solving steps (e.g., for debugging adaptive behavior), you could use a callback function with solve_ivp, but for most use cases (printing/plotting aligned data), the post-solve calculation is simpler and cleaner.

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