%一维谐振子运动方程的数值求解
clear all
close all
clc
m=1; %谐振子质量
k=4*pi*pi;%弹簧劲度系数


h=0.0001; %时间步长,0.5,0.2,0.05
x0=1;   %初始位置,可以任意设置
p0=0;   %初始动量,可以任意设置
a=0;    %时间起点
b=50;   %时间终点
t=a:h:b;
N=(b-a)/h;
x=zeros(1,N+1);
p=zeros(1,N+1);
x(1)=x0;
p(1)=p0;
omega=pi/2;
An=-5000;
for n=1:N %Euler法迭代公式
    x(n+1)=x(n)+h*p(n)/m;
    p(n+1)=p(n)-h*k*x(n)-h*0.6*p(n)+h*An*cos(omega*t(n));
end

subplot(1,2,1),plot(t,x,'r')
title('谐振子的x(t)曲线')
subplot(1,2,2),plot(x,p,'b')
title('谐振子的x-p相图')