相图与相曲线的绘制

这篇文章中，UncleBob 对于使用数学软件进行系统相图与相曲线的绘制方法进行了小结。

相图，相曲线

设$U\subset\mathbb{R}^{n+1}$是开集，$f:U\to\mathbb{R}^n$连续。
对于以下$n$维系统上的IVP（初值问题）

$$\begin{cases} \dot{X}=f(t,X)\\X(t_0)=X_0 \end{cases}$$

它的解$\varphi(t):I\to\mathbb{R}^n$满足

$$\dot{\varphi}(t)=\frac{d\varphi}{dt}(t)=f(t,\varphi(t))$$

且满足初值$\varphi(t_0)=X_0$。

设$\Gamma$是$U$上$C^1$光滑的曲线，具有如下的参数化：

$$\Gamma:=\{(t,\varphi(t)):t\in I\}$$

此时称$\Gamma$是系统$\dot{X}=f(t,X)$的一条积分曲线。

更特殊地，设$V\subset\mathbb{R}^{n}$是开集，$g:V\to\mathbb{R}^d$连续。
我们称系统

$$\dot{X}=g(X)$$

是一个自治系统，$V$为相空间，$g$为向量场（相图）。
设它的解为$\varphi(t):I\to U$，此时称曲线

$$C_\varphi:=\{\varphi(t):t\in I\}$$

是一条相曲线。

单纯看定义可能有些抽象，我们在下面的实际操作观察。

Mathematica实现

系统1

考虑如下系统

$$\dot{X}=JX$$

这里$J=\begin{bmatrix}\lambda&1\\&\lambda\end{bmatrix}$

使用Mathematica绘制相图：

Manipulate[StreamPlot[{Lambda*x + y, Lambda*y}, {x, -5, 5}, {y, -5, 5}], {Lambda, -3, 3}]

这里的Manipulate函数让我们能很方便的交互式操作调整$\lambda$的值观察相曲线的变化。

系统2

考虑如下系统

$$\begin{cases} \dot{x}=\sin x\\\dot{y}=\sin y \end{cases}$$

使用Mathematica绘制相图：

StreamPlot[{Sin[x], Sin[y]}, {x, -5, 5}, {y, -5, 5}]

系统3

考虑如下系统（Lorenz系统）

$$\begin{cases} \dot{x}=\sigma(y-x)\\\dot{y}=rx-y-xz\\\dot{z}=xy-bz \end{cases}$$

取$\sigma=10,b=8/3,r=28$，选取初值$t_0=0,X_0=(11,45,14)$，

使用Mathematica绘制相曲线：

s = NDSolve[{x'[t] == 10 (y[t] - x[t]), y'[t] == 28*x[t] - y[t] - x[t]*z[t], z'[t] == x[t]*y[t] - 8/3*z[t],x[0] == 11, y[0] == 45, z[0] == 14}, {x, y, z}, {t, 50}]
ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. s], {t, 0, 50}, PlotRange -> Full]

Matlab实现

与上述系统相同。

对于第一个：

% 定义参数
lambda = -0.5; % 调整lambda值以查看不同的动态行为
J = [lambda, 1; 0, lambda];

% 定义网格
[x1, x2] = meshgrid(-5:0.5:5, -5:0.5:5);

% 计算方向场
dx1 = lambda*x1 + x2;
dx2 = lambda*x2;

% 绘制方向场
figure;
quiver(x1, x2, dx1, dx2, 'b');
hold on;

title(['Phase Portrait with \lambda = ', num2str(lambda)]);

对于第二个

% 定义网格
[x, y] = meshgrid(-2*pi:0.5:2*pi, -2*pi:0.5:2*pi);

% 定义方向场
dx = sin(x);
dy = sin(y);

% 绘制方向场
figure;
quiver(x, y, dx, dy, 'b');
hold on;

% 绘制一些轨迹
% 定义初始条件
initial_conditions = [-6 -6; -4 -2; 0 0; 2 4; 6 6];
tspan = [0, 10];

for i = 1:size(initial_conditions, 1)
    [t, traj] = ode45(@(t, X) [sin(X(1)); sin(X(2))], tspan, initial_conditions(i, :)');
    plot(traj(:,1), traj(:,2), 'r', 'LineWidth', 1.5);
end

title('Phase Portrait of \dot{x} = sin(x), \dot{y} = sin(y)');

这里额外添加了一些典型相曲线。

对于第三个

 % 参数定义
sigma = 10;
b = 8/3;
r = 28;

% 初始条件和时间范围
X0 = [11, 45, 14]; % 初值
tspan = [0, 50];   % 时间范围

% 定义Lorenz系统
lorenz = @(t, X) [ ...
    sigma * (X(2) - X(1));           % dx/dt
    r * X(1) - X(2) - X(1) * X(3);  % dy/dt
    X(1) * X(2) - b * X(3)          % dz/dt
];

% 使用ode45求解
[t, X] = ode45(lorenz, tspan, X0);

% 绘制相轨迹
figure;
plot3(X(:,1), X(:,2), X(:,3), 'b', 'LineWidth', 1.5);
grid on;

% 图形修饰
xlabel('x');
ylabel('y');
zlabel('z');
title('Lorenz System Phase Curve (\sigma=10, b=8/3, r=28)');
view(3); % 设置三维视角

其他实现

例如对于系统1，使用python中的Matplotlib和NumPy绘制相图

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

# 定义参数
lambda_val = -0.5  # 调整lambda值查看不同行为
J = np.array([[lambda_val, 1], [0, lambda_val]])

# 定义系统的右边函数
def system(t, X):
    return J @ X

# 定义网格
x1 = np.linspace(-5, 5, 20)
x2 = np.linspace(-5, 5, 20)
x1, x2 = np.meshgrid(x1, x2)

# 计算方向场
dx1 = lambda_val * x1 + x2
dx2 = lambda_val * x2

# 归一化方向场向量，可以删除
magnitude = np.sqrt(dx1**2 + dx2**2)
dx1 /= magnitude
dx2 /= magnitude

# 绘制方向场
plt.figure(figsize=(8, 6))
plt.quiver(x1, x2, dx1, dx2, color='b', alpha=0.6)

# 绘制轨迹
initial_conditions = [[-4, -4], [-3, 3], [4, -4], [3, 3]]  # 初始条件
t_span = [0, 10]  # 时间范围

for ic in initial_conditions:
    sol = solve_ivp(system, t_span, ic, t_eval=np.linspace(0, 10, 200))
    plt.plot(sol.y[0], sol.y[1], 'r', linewidth=1.5)

# 图形修饰
plt.title(f'Phase Portrait with λ = {lambda_val}', fontsize=14)
plt.xlabel('x₁', fontsize=12)
plt.ylabel('x₂', fontsize=12)
plt.axhline(0, color='black',linewidth=0.5, linestyle='--')
plt.axvline(0, color='black',linewidth=0.5, linestyle='--')
plt.grid(True)
plt.axis('equal')
plt.show()

这里我们进行了归一化处理（让每个向量长度一致），实际上是否归一化要看不同的场景。