拉格朗日方程的数值解（Matlab）

贡献者： addis

本文处于草稿阶段。

预备知识　哈密顿正则方程，偏导与差分

　　若程序中给出拉格朗日量拉格朗日方程的数值函数 $L (q, \dot{q}, t)$ ，输入和输出均为数值（例如双精度实数），那么如何数值求解运动方程呢？

　　问题的关键是如何列出一个一阶常微分方程组。如果可以给出系统哈密顿量的数值函数，那么数值解哈密顿正则方程是首选的做法，因为它本身已经是一阶常微分方程组，可以直接用解算器进行数值解。但对于较复杂的系统，哈密顿量比拉格朗日量难算得多，甚至可能没有解析表达式。

　　相比之下，拉格朗日量虽然容易写出，但方程的数值解比哈密顿正则方程要难一些，且如果用差分法计算微分会引入一定的数值误差。

　　由拉格朗日方程（式 3 ）

\begin{matrix} (1) & \frac{d}{d t} \frac{\partial L}{\partial {\dot{q}}_{i}} = \frac{\partial L}{\partial q_{i}} + Q_{i} (i = 1, \dots, N) . \end{matrix}

根据全微分，左边有

\begin{matrix} (2) & \frac{d}{d t} \frac{\partial L}{\partial {\dot{q}}_{i}} = \sum_{j} \frac{\partial}{\partial {\dot{q}}_{j}} \frac{\partial L}{\partial {\dot{q}}_{i}} {\ddot{q}}_{j} + \sum_{j} \frac{\partial}{\partial q_{j}} \frac{\partial L}{\partial {\dot{q}}_{i}} {\dot{q}}_{j} + \frac{\partial}{\partial t} \frac{\partial L}{\partial {\dot{q}}_{i}} (i, j = 1, \dots, N) . \end{matrix}

代入得

\begin{matrix} (3) & \sum_{j} \frac{\partial}{\partial {\dot{q}}_{j}} \frac{\partial L}{\partial {\dot{q}}_{i}} {\ddot{q}}_{j} = \frac{\partial L}{\partial q_{i}} - \sum_{j} \frac{\partial}{\partial q_{j}} \frac{\partial L}{\partial {\dot{q}}_{i}} {\dot{q}}_{j} - \frac{\partial}{\partial t} \frac{\partial L}{\partial {\dot{q}}_{i}} + Q_{i} . \end{matrix}

这样，二阶导数只存在于左边的

{\ddot{q}}_{j}

，其他项都只是

q, \dot{q}, t

的函数。该式中的偏微最好计算出解析表达式，但也可以通过差分来计算(会引入更多误差）。数值解线性方程¹，就可以得到（令

v_{i} = {\dot{q}}_{i}

）（定理 2 ）

\begin{matrix} (4) & {\begin{aligned} {\dot{v}}_{i} = f_{i} (q, v_{i}, t) \\ {\dot{q}}_{i} = v_{i} \end{aligned} (i = 1, \dots, N) . \end{matrix}

就得到了

2 N

条方程组成的一阶常微分方程组，变量一共有

2 N

个。可以使用四阶龙格库塔法等方法求解。

　　代码依赖 “偏导与差分” 中的 D2_ij.m，以及 “导数与差分” 的 D_i.m。

未完成：dL, d2L, Q 参数未测试

代码 1：lagr_solve.m

% 数值解拉格朗日方程
% 拉格朗日量 L(qqt), qqt = [q, dq, t]
% 一阶偏微分 dL(i,qqt), 可选, 若不支持 i > N， 返回 nan
% 二阶偏微分 d2L(i,j,qqt) 可选, i = N+1:2*N, j = 1:2*N+1
% 广义力 Q(i, qqt) 可选
function qq = lagr_solve(L, qq0, t, RelTol, dL, d2L, Q, vpa_flag)
%====== 参数设置 =======
h = 1e-5; % 1 阶差分长度
h2 = 1e-5; % 2 阶差分长度
%=====================
N = numel(qq0)/2;
if ~exist('vpa_flag', 'var') || isempty(vpa_flag)
    vpa_flag = false;
end
if ~exist('d2L', 'var') || isempty(d2L)
    if ~exist('dL', 'var') || isempty(dL) || isnan(dL(N+1,[qq0, tmin]))
        if vpa_flag
            d2L = @(i,j,qqt)D2_ij_vpa(L,i,j,qqt,h2);
        else
            d2L = @(i,j,qqt)D2_ij(L,i,j,qqt,h2);
        end
    else
        d2L = @(i,j,qqt)d2L_dif1(dL,i,j,qqt,h);
    end
end
if ~exist('dL', 'var') || isempty(dL)
    if vpa_flag
        dL = @(i, qqt)D_i_vpa(L, i, qqt, h);
    else
        dL = @(i, qqt)D_i(L, i, qqt, h);
    end
end
if ~exist('Q', 'var') || isempty(Q)
    Q = @(i,qqt) 0;
end
options = odeset('RelTol', RelTol);
[T, QQ] = ode45(@(t,qq)ode_fun(t, qq, dL, d2L, Q), ...
        [t(1), t(end)], qq0, options);
qq = zeros(numel(t), 2*N);
for i = 1:2*N
    qq(:,i) = interp1(T, QQ(:,i), t, 'spline');
end
end

function dqq = ode_fun(t, qq, dL, d2L, Q)
N = numel(qq)/2;
dqq = zeros(2*N, 1);
dqq(1:N) = qq(N+1:end);
qqt = [qq; t];
A = zeros(N, N); % A_ij = d2L/d(\dot q_i)d(\dot q_j)
for i = 1:N
    for j = i:N
        A(i,j) = d2L(i+N, j+N, qqt);
        if i ~= j
            A(j,i) = A(i,j);
        end
    end
end
b = zeros(N, 1);
for i = 1:N
    b(i) = dL(i, qqt) - d2L(i+N,2*N+1,qqt) + Q(i,qqt);
    for j = 1:N
        b(i) = b(i) - d2L(N+i,j,qqt)*qqt(j+N);
    end
end
dqq(N+1:end) = A \ b;
if any(isnan(dqq) | isinf(dqq))
    error('wrong!');
end
end

% d^2 L / d(qqt_i)d(qqt_j)
% dL is dL/d(qqt_i)
function ret = d2L_dif1(dL, i, j, qqt, h)
qqt(j) = qqt(j) - 0.5*h;
L1 = dL(i, qqt, h);
h = (qqt(j) + h) - qqt(j);
qqt(j) = qqt(j) + h;
L2 = dL(i, qqt, h);
ret = (L2 - L1) / h;
end

　　用双摆的例子来测试，比较 “双摆的数值计算（Matlab）” 的结果：

代码 2：lagr_solve_demo.m

% lagr_solve_demo
function lagr_solve_demo
% ==== 参数 =====
m1 = 1; m2 = 1; g = 9.8;
l1 = 1; l2 = 1;
tmin = 0; tmax = 10; Nt = 200;
qq0 = [pi, pi, -4, 6];
RelTol = 1e-6;
L = @(qqt)L_fun(qqt, m1, l1, m2, l2, g);
% ==============
close all;
t = linspace(tmin, tmax, Nt);
qq = lagr_solve(L, qq0, t, RelTol);

th1 = qq(:,1); th2 = qq(:,2);
x1 = l1*sin(th1); y1 = -l1*cos(th1);
x2 = x1 + l2*sin(th2); y2 = y1 - l2*cos(th2);
figure;
% 动画
for it = 1:Nt
    clf;
    plot([0, x1(it)], [0, y1(it)], 'k'); hold on;
    scatter(x1(it), y1(it), 'k');
    plot([x1(it), x2(it)], [y1(it), y2(it)], 'k');
    scatter(x2(it), y2(it), 'k');
    plot(x2(1:it), y2(1:it));
    axis equal; axis([-1,1,-1,1]*(l1+l2)*1.1);
    title(['t = ' num2str(t(it), '%.2f')]);
    xlabel x; ylabel y;
    pause(0.1); % saveas(gcf, [num2str(it) '.png']);
end
end

% 双摆
function ret = L_fun(qqt, m1, l1, m2, l2, g)
th1 = qqt(1); th2 = qqt(2);
w1 = qqt(3); w2 = qqt(4);
T = 0.5*m1*(l1*w1)^2 + 0.5*m2*((l1*w1*cos(th1) + l2*w2*cos(th2))^2 ...
    + (l1*w1*sin(th1) + l2*w2*sin(th2))^2);
V = -m1*g*l1*cos(th1) - m2*g*(l1*cos(th1) + l2*cos(th2));
ret = T - V;
end

1. ^ 如果 $N$ 不大可以直接用克拉默法则，甚至写出解析解。

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。