假设有一恒星系这个恒星系统中有一颗行星而这个行星有一颗卫星假设恒星重量占恒星系总质量的百分之9999因此可以认为恒星不动。用拉格朗日方程求解并用matlab画出运动轨迹

由于题目中只给出了恒星的质量占比，没有具体数值，因此我们无法得出行星和卫星的质量，也无法画出具体的运动轨迹。下面给出一般情况下的求解方法和matlab代码。

假设行星质量为$m_1$，卫星质量为$m_2$，行星和卫星之间的距离为$r$，行星和恒星之间的距离为$d$，卫星和恒星之间的距离为$l$。建立以恒星为原点的坐标系，选取行星和卫星之间连线的中点为原点。则行星和卫星分别沿着径向和切向运动，其运动轨迹可以用极坐标表示。设行星的极角为$\theta_1$，卫星的极角为$\theta_2$，则有：

$$ r = \frac{d + l}{2}, \quad \theta_1 = \theta_2 - \frac{\pi}{2} $$

根据拉格朗日方程，有：

$$ \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}}\right) - \frac{\partial L}{\partial \theta} = 0 $$

其中$L$为系统的拉格朗日量，可以表示为：

$$ L = \frac{1}{2}m_1(\dot{r}^2 + r^2\dot{\theta_1}^2) + \frac{1}{2}m_2(\dot{r}^2 + r^2\dot{\theta_2}^2) - \frac{Gm_1m_2}{d} - \frac{Gm_1M}{l} - \frac{Gm_2M}{l} $$

其中$G$为万有引力常数，$M$为恒星质量。代入$r$和$\theta_1$的表达式，可以得到：

$$ L = \frac{1}{2}m_1\left(\dot{d}^2 + \frac{(d+l)^2}{4}\dot{\theta_2}^2\right) + \frac{1}{2}m_2\left(\dot{d}^2 + \frac{(d+l)^2}{4}\dot{\theta_2}^2\right) - \frac{Gm_1m_2}{d} - \frac{Gm_1M}{\sqrt{(d+l)^2/4+l^2}} - \frac{Gm_2M}{\sqrt{(d+l)^2/4+l^2}} $$

化简可得：

$$ L = \frac{1}{2}m_1\dot{d}^2 + \frac{1}{8}(m_1+m_2)(d+l)^2\dot{\theta_2}^2 - \frac{Gm_1m_2}{d} - \frac{Gm_1M}{\sqrt{(d+l)^2/4+l^2}} - \frac{Gm_2M}{\sqrt{(d+l)^2/4+l^2}} $$

代入拉格朗日方程，可得：

$$ \ddot{d} = \frac{G(m_1+m_2)}{d^2} - \frac{Gm_1}{(d+l)^2/4+l^2)^{3/2}} - \frac{Gm_2}{(d+l)^2/4+l^2)^{3/2}} $$

$$ \ddot{\theta_2} = -\frac{2\dot{d}\dot{\theta_2}}{d+l} - \frac{Gm_1}{(d+l)^2/4+l^2)^{3/2}}\sin\theta_2 - \frac{Gm_2}{(d+l)^2/4+l^2)^{3/2}}\sin(\theta_2 - \theta_1) $$

这是一个二阶常微分方程组，可以用matlab的ode45函数求解。先定义一个函数，用于返回方程组的右端项：

function ddt = orbit(t, y, m1, m2, G, l)
    d = y(1);
    theta2 = y(2);
    dd = y(3);
    dtheta2 = y(4);
    
    r = (d + l) / 2;
    theta1 = theta2 - pi/2;
    
    ddd = G*(m1+m2)/d^2 - G*m1/((d+l)^2/4+l^2)^1.5 - G*m2/((d+l)^2/4+l^2)^1.5;
    dtheta2d = -2*dd*dtheta2/(d+l) - G*m1/((d+l)^2/4+l^2)^1.5*sin(theta2) - G*m2/((d+l)^2/4+l^2)^1.5*sin(theta2-theta1);
    
    ddt = [dd; dtheta2; ddd; dtheta2d];
end

然后设定初始条件和参数，调用ode45函数求解：

% 初始条件和参数
d0 = 1;
theta20 = 0;
dd0 = 0;
dtheta20 = 1;

m1 = 1;
m2 = 0.1;
M = 100;
G = 1;

l = 2;

tspan = [0, 100];
y0 = [d0; theta20; dd0; dtheta20];

% 求解方程
[t, Y] = ode45(@(t, y)orbit(t, y, m1, m2, G, l), tspan, y0);

% 画图
figure
polarplot(Y(:,2), Y(:,1))
``