含依赖参数的ODE求解与相位平面/箭头图绘制技术求助

Hey there! I see you're working on solving a linear ODE system and generating phase portraits for different k values. Your current code is a great start, but we need to tweak it to handle multiple k's cleanly and add the vector field (arrow plot) for each phase plane. Let's break this down step by step:

1. Fix the Model Function to Accept k as a Parameter #

Right now, your model function relies on a global k, which gets messy when looping through multiple values. We'll modify it to take k as an explicit argument, which we can pass to odeint using the args parameter.

2. Generate Vector Fields for Phase Portraits #

To create the arrow plot (vector field), we'll set up a grid of (x, y) points, calculate the direction of the ODE at each point, and plot them with plt.quiver. We'll normalize the vectors so arrow sizes stay consistent across the plot.

3. Loop Through Each k Value #

We'll define a list of k values to test, then for each k:

• Solve the ODE using your initial condition
• Plot the trajectory on the phase plane
• Overlay the vector field
• Add labels and a k-specific title for clarity

Full Working Code #

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

def model(X, t, k):
    x, y = X
    dxdt = k * x - y
    dydt = x + y
    return [dxdt, dydt]

# Initial state
X0 = [1, 1]

# Time range for integration
t = np.linspace(0, 10, 1000)

# List of k values to test (adjust this as needed)
k_values = [-2, 0, 2]

# Set up subplots: one for each k value
fig, axes = plt.subplots(nrows=len(k_values), ncols=1, figsize=(6, 4*len(k_values)))
fig.tight_layout(pad=3.0)

# Loop through each k and generate its phase portrait
for idx, k in enumerate(k_values):
    ax = axes[idx] if len(k_values) > 1 else axes
    
    # Solve ODE for the current k
    X = odeint(model, X0, t, args=(k,))
    x_trajectory = X[:, 0]
    y_trajectory = X[:, 1]
    
    # Create grid for vector field
    x_grid, y_grid = np.meshgrid(np.linspace(-5, 5, 20), np.linspace(-5, 5, 20))
    dx_grid = k * x_grid - y_grid
    dy_grid = x_grid + y_grid
    
    # Normalize vectors to keep arrow size consistent
    magnitude = np.sqrt(dx_grid**2 + dy_grid**2)
    dx_normalized = dx_grid / magnitude
    dy_normalized = dy_grid / magnitude
    
    # Plot the vector field (light gray arrows)
    ax.quiver(x_grid, y_grid, dx_normalized, dy_normalized, color='lightgray', alpha=0.6)
    
    # Plot the trajectory from initial condition
    ax.plot(x_trajectory, y_trajectory, color='darkblue', linewidth=2, label=f'Trajectory (X0={X0})')
    
    # Add plot details
    ax.set_title(f'Phase Portrait for k = {k}')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-5, 5)
    ax.set_ylim(-5, 5)

plt.show()

Key Notes: #

• Parameterized Model: The model function now accepts k as an argument, making it easy to test different values without redefining the function.
• Vector Field: We use meshgrid to create a grid of points, compute the ODE's direction at each point, and normalize vectors so arrows don't grow too large in high-magnitude regions.
• Subplots: Each k gets its own subplot, arranged vertically. Adjust nrows/ncols if you prefer a horizontal layout.
• Trajectory: The solution from odeint is plotted as a dark blue line over the vector field, showing how the system evolves from your initial condition [1,1].

Feel free to tweak the k_values list, time range, grid limits, or initial condition to match your specific needs!

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