在matlab中通过哈密顿正则方程的零本征值以及非零本征值对应的本征向量组成的解来做的数值算例代码

以下是一个使用MATLAB求解哈密顿正则方程的数值算例代码：

% 定义哈密顿量
H = @(p,q) 0.5*p^2 + 0.5*q^2 + 0.1*q^4;

% 定义哈密顿正则方程
dqdt = @(p) p;
dpdt = @(q) -q - 0.4*q^3;

% 定义时间步长和总时间
dt = 0.01;
T = 20;

% 定义初始条件
q0 = 0.5;
p0 = 1;

% 初始化解向量
t = 0:dt:T;
q = zeros(size(t));
p = zeros(size(t));
q(1) = q0;
p(1) = p0;

% 使用龙格库塔方法求解哈密顿正则方程
for i = 2:length(t)
    k1q = dqdt(p(i-1));
    k1p = dpdt(q(i-1));
    k2q = dqdt(p(i-1) + 0.5*dt*k1p);
    k2p = dpdt(q(i-1) + 0.5*dt*k1q);
    k3q = dqdt(p(i-1) + 0.5*dt*k2p);
    k3p = dpdt(q(i-1) + 0.5*dt*k2q);
    k4q = dqdt(p(i-1) + dt*k3p);
    k4p = dpdt(q(i-1) + dt*k3q);
    q(i) = q(i-1) + (1/6)*dt*(k1q + 2*k2q + 2*k3q + k4q);
    p(i) = p(i-1) + (1/6)*dt*(k1p + 2*k2p + 2*k3p + k4p);
end

% 绘制轨迹
figure;
plot(q,p);
xlabel('q');
ylabel('p');
title('Trajectory');

% 计算本征值和本征向量
[V,D] = eig([0 1; -1 -0.4*q(end)^2]);

% 绘制本征向量
hold on;
quiver(0,0,V(1,1),V(2,1),'LineWidth',2);
quiver(0,0,V(1,2),V(2,2),'LineWidth',2);
legend('Trajectory','Eigenvalue 1','Eigenvalue 2');

这个代码求解了一个简单的二维哈密顿系统，并使用龙格库塔方法求解哈密顿正则方程，最后计算了系统的本征值和本征向量，并将它们绘制在轨迹上。注意，由于哈密顿系统的能量守恒性质，轨迹应该在哈密顿量等于初始条件的值的等能线上运动。