%例题2：一维谐振子运动方程的数值求解
clear all
close all
clc
m=1; %谐振子质量
k=0.1;%弹簧劲度系数


h=0.001; %时间步长,0.5,0.2,0.05
x0=5;   %初始位置,可以任意设置
p0=0;   %初始动量,可以任意设置
a=0;    %时间起点
b=200;   %时间终点
t=a:h:b;
N=(b-a)/h;
x=zeros(1,N+1);
p=zeros(1,N+1);
x(1)=x0;
p(1)=p0;

for n=1:N %Euler法迭代公式
    x(n+1)=x(n)+h*p(n)/m;
    p(n+1)=p(n)-h*k*x(n);
end
%%%%%%%%%%解析结果%%%%%%%%%%%%%%%%%%%%
A=sqrt(x0^2+p0^2/(m*k)); %谐振子的振幅
pha=acos(x0/A);          %谐振子的初始相位，若x0与p0同号，pha为负，否则pha为正
xx=A*cos(sqrt(k/m)*t+pha); %谐振子位置的精确解
pp=-m*A*sqrt(k/m)*sin(sqrt(k/m)*t+pha); %动量的精确解
EE=pp.^2/(2*m)+0.5*k*xx.^2; %能量的精确值

%%%%%%%%%%Euler法数值结果%%%%%%%%%%%%%%%%%%%%
E=p.^2/(2*m)+0.5*k*x.^2; %能量的数值解
subplot(3,1,1),plot(t,x,'r',t,xx,'b')
title('谐振子的x(t)曲线')
legend('Euler法','精确解')
subplot(3,1,2),plot(x,p,'r',xx,pp,'b')
title('谐振子的x-p相图')
subplot(3,1,3),plot(t,E,'r',t,EE,'b')
title('谐振子的能量随时间变化关系')