众所周知，行星轨道是个封闭的以太阳为焦点的椭圆。这里库伦势场的一个特点是，轨道封闭。按照bertrand定理，轨道封闭的中心势场其实只有两种，即库伦势和谐振子势。

我们可以用四阶runnge-kutta法做下模拟。

考虑中心势，V(r)=-C/r^alpha。其中alpha是关键参数。alpha=1，即库伦势。alpha=-2，即谐振子势。

RK方法不是专门为哈密顿系统设计的算法，所以它没充分考虑哈密顿系统的特性。具体而言，它不保持系统辛结构。各位同学可查看下，在四阶Runge-Kutta法下，系统能量是否保持不变。

RK方法是19世纪计算数学的巅峰之作，考试必考！当然，在电子计算机发明前，计算数学没有充分发展，但即便如此，RK方法依然毫不逊色于20世纪那十大算法。

程序如下：

% 4th-RK method; -C/r^alpha potential

clear all; close all; clc;tic;

C = 1;

alpha = 1.00;

N = 2e4;

ylist = zeros(4, N);

Elist = zeros(1, N);

x0 = [1;0];             % initial position

v0 = [0.5;0.5];         % initial velocity

y = [x0; v0];

ylist(:,1) = y;

Elist(1,1) = 0.5*norm(v0)^2 - C/norm(x0)^alpha;

dt = 0.001;

h1= plot(nan);

hold on

plot(0,0,'O')

plot(x0(1), x0(2),'*')

xlim([-1.5 1.5])

ylim([-1.5 1.5])

xlabel('x','fontsize',20)

ylabel('y','fontsize',20)

str = strcat('alpha=',num2str(alpha));

title(str,'fontsize',20)

axis square

for s = 2: N

   ytemp = y;

   r = sqrt(ytemp(1)^2 + ytemp(2)^2);

   K1 = dt*[ytemp(3); ytemp(4); -alpha*C*ytemp(1)/r^(alpha+2); -alpha*C*ytemp(2)/r^(alpha+2)];

   ytemp = y + 0.5*K1;

   r = sqrt(ytemp(1)^2 + ytemp(2)^2);

   K2 = dt*[ytemp(3); ytemp(4); -alpha*C*ytemp(1)/r^(alpha+2); -alpha*C*ytemp(2)/r^(alpha+2)];

   ytemp = y + 0.5*K2;

   r = sqrt(ytemp(1)^2 + ytemp(2)^2);

   K3 = dt*[ytemp(3); ytemp(4); -alpha*C*ytemp(1)/r^(alpha+2); -alpha*C*ytemp(2)/r^(alpha+2)];

   ytemp = y + K3;

   r = sqrt(ytemp(1)^2 + ytemp(2)^2);

   K4 = dt*[ytemp(3); ytemp(4); -alpha*C*ytemp(1)/r^(alpha+2); -alpha*C*ytemp(2)/r^(alpha+2)];

   y = y + (K1+ 2*K2 + 2*K3 + K4)/6;

   ylist(:,s) = y;

   if mod(s,20)==0

       set(h1,'XData',ylist(1,1:s),'YData',ylist(2,1:s))

       drawnow

       pause(0.001)

   end

   Elist(1,s) = 0.5*norm(y(3:4))^2 - C/norm(y(1:2))^alpha;

end

toc

h2 = figure;

plot(dt*(0:N-1),Elist)

不同alpha值，轨道的开闭情况：