MATLAB 繪製二階微分方程組的 x-y 相軌跡：從數值解到相平面曲線

問題數學模型

給定二階常微分方程組：

y'' = -10.01476 + 0.20811 * (0.0587 * x' - 0.9983 * y')
x'' = -0.20811 * (0.0587 * y' + 0.9983 * x')

初始條件：t=0, x=0, y=10600。目標：直接繪製 x-y 相平面曲線，而非 x-t、y-t 時間序列。

解題思路：降階 + 數值積分

1. 引入速度變數：vx = x', vy = y'，將二階系統轉為四階一階系統：

dx/dt  = vx
dy/dt  = vy
dvx/dt = -0.20811 * (0.0587*vy + 0.9983*vx)
dvy/dt = -10.01476 + 0.20811 * (0.0587*vx - 0.9983*vy)

2. 使用 ode45 求解：得到 t, x(t), y(t), vx(t), vy(t)。
3. 繪製相軌跡：plot(x, y) 即為所求 x-y 曲線。

MATLAB 完整程式碼

function dxdt = odefun(t, X)
    x  = X(1); y  = X(2);
    vx = X(3); vy = X(4);
    dxdt = zeros(4,1);
    dxdt(1) = vx;
    dxdt(2) = vy;
    dxdt(3) = -0.20811 * (0.0587*vy + 0.9983*vx);
    dxdt(4) = -10.01476 + 0.20811 * (0.0587*vx - 0.9983*vy);
end

% 主程式
X0 = [0; 10600; 0; 0];          % 假設初速為 0，若非請修正
tspan = [0 100];                % 依物理意義調整積分時間
[t, X] = ode45(@odefun, tspan, X0);

x = X(:,1); y = X(:,2);
figure; plot(x, y, 'b-', 'LineWidth', 1.5);
xlabel('x'); ylabel('y'); title('x-y 相平面軌跡'); grid on;

進階：直接求 dy/dx 關係 (解析嘗試)

若不需時間參數，可嘗試消去 t：

dy/dx = vy/vx
d²y/dx² = (vx*dvy/dt - vy*dvx/dt) / vx^3

但非線性耦合強，解析解極難，**數值法為首選**。

常見坑點

• 初始速度未給：題目僅給位置，若 vx(0), vy(0) 不為零，曲線形狀將劇烈不同，需確認物理語境。
• 積分區間過短：tspan 需涵蓋系統穩態或週期特徵，建議先大區間粗看再縮小。
• 硬式方程：若參數導致剛性，改用 ode15s 或 ode23s。

結果驗證

執行上述代碼將直接輸出 x-y 相平面曲線**，滿足「不要 x-t、y-t，直接要 x-y 曲線」的需求。若需匯出座標點，使用 writematrix([x y], 'trajectory.csv') 即可。