本轮均轮背后的傅里叶分解原理（matlab演示）

1.简介

均轮、本轮是希腊天文学家希帕克斯提出的以地球为中心的学说，认为天体沿着本轮做匀速圆周运动，这些本轮的中心又沿着各自的更大的均轮做以地球为中心的匀速圆周运动。

在长达十四个世纪的古代时期，占据主导地位的天文学理论，是托勒密的地心说。由于地心说中并不是所有天体都按照围绕地球，做圆形选转的运动方式来运动，所以作为修正，提出来本轮和均轮的概念。

比如上图中的行星视运动中“顺-留-逆-留-顺”运行轨迹，是本轮运动和均轮运动两者综合的结果。

这个方法可以十分精准的预测行星的位置，但是代价便是越来越多的新的本轮均轮被加原来均轮上，形成大圆套小圆的结构。甚至有的天体用到了12个圆来进行模拟。

在当然根据现代的观测结果，人们发现将坐标系放在太阳上可以有效的简化太阳系内行星的运动轨迹，或者说把所有因太阳的引力而绕太阳为焦点的所有行星称为太阳系的行星。此外根据坐标变换原则，人们还可以构建各种奇葩新颖的坐标系，但这不是本文的讨论重点。

本文利用傅里叶分解的原理，证明所有的平面运动都可以用本轮-均轮系统表示，而且圆越多精度越高。而且也存在有限个圆就能完全拟合的曲线。

2.用到的数学工具（傅里叶分解）

这里用到了傅里叶分解，由于研究的是二维曲线，所以需要进行x和y的两次分解。然而这里我直接用复数形式表示，这样求解一次就可以直接出现xy的方程。而且附带的好处是一次傅里叶分解如果出现角速度大小相等方向相反的时候，会自动合并，而复数域分解就会自动分成正负两项，提供了后续绘图的便利。

复数域的傅里叶分解公式：
$$f(x)=\sum _{n=-\infty }^{\infty }{c}_n\cdot exp\frac{in\pi x}{l}$$
$$c_n=\frac{1}{2l}{\int }_l^{-l}f(t)\cdot exp\frac{-in\pi t}{l}dt$$

3.绘图演示

3.1 椭圆

对于长轴Lx=2,短轴Ly=1的椭圆，可以傅里叶分解后出现2项，且其余高阶项均为0.也就是说，椭圆用2个圆就可以完全精确的拟合出来。
圆图：

棍图：

3.2 三角形

三角形由于存在尖角，所以不能用有限多的圆形去精确拟合，但是可以增加圆形来逐渐接近。

例：2个圆形的拟合

例：4个圆形的拟合

随着圆的拟合越多，可以看到拟合越来越接近。而且在本例题的实际操作中发现偶数个圆比奇数个圆拟合效果更好。

3.3 任意图形

这里为了突出本轮均轮系统的强悍性，任意的图形都可以进行分解。

其中要求为x区间是沿y轴对称的。f虽然为x的参数，但是不一定要写出关系式，只要是一 一对应的就行。而且f最好是收尾封闭的，虽然不封闭也能求解出结果，但效果不好。f平面图形曲率光滑，拟合效果一般会更好。f图形曲线可以出现交叉。
f的起点最好取在x轴上，f的型心最好在坐标轴原点上，在计算过程中可以减少某些项，但是由于这是编程实现的，所以这个要求也不是必要的。

我在这里随便画了一个二维曲线，傅里叶分解如下：

下图：实际s=5的理想曲线（红色）与拟合曲线（蓝色）↓

下图：圆图↓

下图：棍图↓

4.matlab代码

clear

%要求，1给出精度dx
%2对称单增参数x
%3复数域函数f，收尾相接，不一定与x有特别的关系
%4给出拟合圆的阶数s

dx=1e-3;
%例1，三角形
% x1=-3/2:dx:-1/2-dx;
% x2=-1/2:dx:1/2-dx;
% x3=1/2:dx:3/2;
% x=[x1,x2,x3];
% f1=(-5/2-3*x1)+1i*(3*sqrt(3)/2+sqrt(3)*x1);
% f2=(-1)+1i*(-2*sqrt(3)*x2);
% f3=(-5/2+3*x3)+1i*(-3*sqrt(3)/2+sqrt(3)*x3);
% f=[f1,f2,f3];

% %例2，椭圆
% x=0:dx:2*pi;
% f=2*cos(x)+1i*sin(x);

%例3，多边形
Nx=1/dx;
f=[linspace(0,5,Nx)+1i*linspace(0,0,Nx),...
    linspace(5,5,Nx)+1i*linspace(0,1,Nx),...
    linspace(5,3,Nx)+1i*linspace(1,1,Nx),...
    linspace(3,3,Nx)+1i*linspace(1,5,Nx),...
    linspace(3,0,Nx)+1i*linspace(5,2,Nx),...
    linspace(0,2,Nx)+1i*linspace(2,2,Nx),...
    linspace(2,2,Nx)+1i*linspace(2,1,Nx),...
    linspace(2,0,Nx)+1i*linspace(1,1,Nx),...
    linspace(0,0,Nx)+1i*linspace(1,0,Nx)];
x=0:length(f)-1;
x=x/max(x)*5;
x=x-max(x)/2;



L=(x(end)-x(1))/2;
s=5;%拟合圆的阶数
N=[-s:s];
%求解半径Cn
Cn=zeros(length(N),1);
for j=1:length(N)
    n=N(j);
    Cn(j)=1/2/L*(trapz(x,f.*exp(-1i*n*pi*x/L)));
end
%化简，根据dx把Cn的小量删掉
Cn_r=real(Cn);Cn_r(abs(Cn_r)<dx)=0;
Cn_i=imag(Cn);Cn_i(abs(Cn_i)<dx)=0;%一般这个都是0
Cn=Cn_r+1i*Cn_i;


%求解轨迹Zn
Zn=zeros(length(N),length(x));%Zn的格式为0，1正负，2正负，3正负。。。
for j=0:s
    if j==0
        Zn(j+1,:)=Cn(s+1)*exp(1i*j*pi*x/L);
    else
        Zn(j*2,:)=Cn(s+1+j)*exp(1i*j*pi*x/L);%角速度为正
        Zn(j*2+1,:)=Cn(s+1-j)*exp(-1i*j*pi*x/L);%角速度为负
    end
end
Znsum=cumsum(Zn,1);


%绘制动图
figure(1)
pic_num=1;
for j=1:40:length(x)
    clf
    hold on
    %绘制圆
    for k=1:s
        %圆心 Zn(k-1)
        %半径 Cn
        %先画正的
        if Cn(s+1+k)~=0
            pos=[real(Znsum(2*k-1,j))-abs(Cn(s+1+k)),imag(Znsum(2*k-1,j))-abs(Cn(s+1+k)),abs(2*Cn(s+1+k)),abs(2*Cn(s+1+k))];
            rectangle('Position',pos,'Curvature',[1 1])
        end
        %再画负的
        if Cn(s+1-k)~=0
            pos=[real(Znsum(2*k,j))-abs(Cn(s+1-k)),imag(Znsum(2*k,j))-abs(Cn(s+1-k)),abs(2*Cn(s+1-k)),abs(2*Cn(s+1-k))];
            rectangle('Position',pos,'Curvature',[1 1])
        end
        
    end
    %绘制曲线
    plot(Znsum(end,1:j),'-r')
    hold off
    %xlim([-3,3]);ylim([-3,3])
    xlim([-1,6]);ylim([-1,6])%画框大小
    pause(0.01)

end

%绘制静图
figure(2)
plot(sum(Zn,1));
hold on
plot(f);
hold off

%绘制棍线动图
pause(0.1)
figure(3)

for j=1:40:length(x)
    clf
    hold on
    %绘制直线
    h=1;
    for k=1:s
        %圆心 Zn(k-1)
        %下一点 Zn(k)
        %先画正的
        if Cn(s+1+k)~=0
            plot([real(Znsum(2*k-1,j)),real(Znsum(2*k,j))],[imag(Znsum(2*k-1,j)),imag(Znsum(2*k,j))],'k-')%,'color',mycolor(colorslect(h),:)
            h=h+1;
        end
        %再画负的
        if Cn(s+1-k)~=0
            plot([real(Znsum(2*k,j)),real(Znsum(2*k+1,j))],[imag(Znsum(2*k,j)),imag(Znsum(2*k+1,j))],'k-')%,'color',mycolor(colorslect(h),:)
            h=h+1;
        end
        
    end
    %绘制曲线
    plot(Znsum(end,1:j),'-r')
    hold off
    %xlim([-2,2]);ylim([-2,2])
    xlim([-1,6]);ylim([-1,6])%画框大小
    pause(0.01)
    

end