洛伦兹方程：100 阶泰勒截断法

下面采用 MATLAB + Symbolic Math Toolbox 实现。将代码保存为 lorenz_taylor_symbolic.m，在 MATLAB R2016b 或更新版本中运行。

参数约定：M = 100 表示泰勒截断阶数；K = 200 表示十进制工作精度，不是时间步数。

采用标准洛伦兹参数 σ = 10、ρ = 28、β = 8/3，步长 h = 1/100 = 0.01。

初始条件：(x₀, y₀, z₀) = (1 + δ, 1, 1)。你没有指定偏置大小，下面默认 δ = 10⁻⁶⁰；未指定终止时间，默认 T = 1（100 步）。两者均可修改。

使用 sym 构造精确常量和初值，通过 vpa 进行 200 位运算。脚本同时计算无扰动参考轨道，方便观察微扰增长。

方程与泰勒系数递推公式

x' = σ(y − x)
y' = ρx − y − xz
z' = xy − βz

每一步在当前时刻 t 展开：
x(t+s) ≈ Σ[n=0,…,M] a_n s^n
y(t+s) ≈ Σ[n=0,…,M] b_n s^n
z(t+s) ≈ Σ[n=0,…,M] c_n s^n

其中 a_n = x^(n)(t)/n!，b_n、c_n 类似。

a_(n+1) = σ(b_n − a_n)/(n+1)

b_(n+1) = [ρa_n − b_n − Σ[j=0,…,n] a_j c_(n−j)]/(n+1)

c_(n+1) = [Σ[j=0,…,n] a_j b_(n−j) − βc_n]/(n+1)

n = 0,…,M−1；最后取 s = h，用 Horner 法计算下一步。

系数递推避免了直接展开庞大的 100 阶符号导数。这里是“符号对象上的高精度求值”，不是将整条轨道保留为纯符号表达式。

%% lorenz_taylor_symbolic.m
% Lorenz系统：M阶泰勒截断 + K位可变精度运算
% 需要 Symbolic Math Toolbox。
% U：无扰动参考轨道；V：仅x初值有微扰的轨道。
clear; clc;

%% 1. 设置参数
M = 100;                   % 泰勒截断阶数
K = 200;                   % 十进制工作精度
digits(K);

% 精确构造参数：避免先经过double近似
sigma = sym(10);
rho   = sym(28);
beta  = sym(8) / sym(3);
h     = sym(1) / sym(100);  % 精确的0.01

delta = sym(10)^(-60);      % 默认微扰，可修改
Tend  = sym(1);             % 默认积分到t=1，可修改

nExact = Tend / h;
assert(isAlways(nExact > 0) && ...
       isAlways(nExact == floor(nExact)), ...
       'Tend必须是h的正整数倍。');

N = double(nExact);        % 仅循环次数转double，不参与状态运算
t = sym((0:N).') * h;       % 精确时间网格

%% 2. 精确初值
u0Exact = [sym(1); sym(1); sym(1)];
v0Exact = [sym(1) + delta; sym(1); sym(1)];

% 不要写 sym(1 + 1e-60)！
% 因为括号内先按double计算，微扰会在转sym之前丢失。
u = vpa(u0Exact, K);
v = vpa(v0Exact, K);

assert(~isAlways(v(1) == u(1)), ...
       '微扰在当前工作精度下已丢失，请提高K或增大delta。');

U = sym(zeros(N+1, 3));
V = sym(zeros(N+1, 3));
U(1,:) = u.';
V(1,:) = v.';

fprintf('M=%d, K=%d, 总步数=%d\n', M, K, N);
fprintf('实际保存的初始x差值：%s\n', ...
        char(vpa(v(1)-u(1), 40)));

%% 3. 时间推进
% 两条轨道分别推进；积分过程不转换为double。
reportEvery = max(1, floor(N/10));
tic;

for k = 1:N
    u = lorenzTaylorStep(u, M, K, h, sigma, rho, beta);
    v = lorenzTaylorStep(v, M, K, h, sigma, rho, beta);

    U(k+1,:) = u.';
    V(k+1,:) = v.';

    if mod(k, reportEvery) == 0 || k == N
        fprintf('已完成 %d/%d 步\n', k, N);
    end
end

elapsed = toc;
fprintf('计算耗时：%.2f 秒\n', elapsed);

%% 4. 高精度差值与最终结果
% 必须先用高精度相减，不能先double(V)-double(U)。
D = vpa(V - U, K);
err = vpa(sqrt(sum(D.^2, 2)), K);
logErr = vpa(log(err) / log(sym(10)), K);

fprintf('\n扰动轨道在 T=%s 时的结果：\n', char(Tend));
names = {'x', 'y', 'z'};
for j = 1:3
    fprintf('%s(T) = %s\n', names{j}, char(vpa(V(end,j), K)));
end
fprintf('\n最终两条轨道的距离：%s\n', char(vpa(err(end), 50)));

%% 5. 绘图：仅显示阶段转double
Ud = double(U);
Vd = double(V);
td = double(t);

figure('Color', 'w');

subplot(1,2,1);
plot3(Ud(:,1), Ud(:,2), Ud(:,3), 'b--');
hold on;
plot3(Vd(:,1), Vd(:,2), Vd(:,3), 'r-');
grid on;
xlabel('x'); ylabel('y'); zlabel('z');
title('Lorenz trajectories');
legend('Unperturbed', 'Perturbed', 'Location', 'best');

subplot(1,2,2);
plot(td, double(logErr), 'k-', 'LineWidth', 1.4);
grid on;
xlabel('t');
ylabel('log_{10} ||u_{\delta}(t)-u_0(t)||_2');
title('Perturbation growth');

%% 6. 保存符号对象，保留高精度结果
save('lorenz_taylor_result.mat', ...
     'M', 'K', 'sigma', 'rho', 'beta', ...
     'h', 'delta', 'Tend', 't', ...
     'u0Exact', 'v0Exact', 'U', 'V', 'D', 'err');

%% 局部函数：单步M阶泰勒法
function next = lorenzTaylorStep(u, M, K, h, sigma, rho, beta)
    % a(n+1)对应a_n = x^(n)(t)/n!，b、c同理。
    % 注意：系数已经包含1/n!，最终求和时不要再除阶乘。
    a = sym(zeros(M+1, 1));
    b = sym(zeros(M+1, 1));
    c = sym(zeros(M+1, 1));

    a(1) = u(1);
    b(1) = u(2);
    c(1) = u(3);

    for n = 0:M-1
        % 非线性乘积的幂级数卷积系数
        % xz_n = sum_{j=0}^n a_j*c_(n-j)
        % xy_n = sum_{j=0}^n a_j*b_(n-j)
        xz_n = vpa(sum(a(1:n+1) .* flipud(c(1:n+1))), K);
        xy_n = vpa(sum(a(1:n+1) .* flipud(b(1:n+1))), K);

        den = sym(n+1);

        a(n+2) = vpa( ...
            sigma*(b(n+1)-a(n+1))/den, K);

        b(n+2) = vpa( ...
            (rho*a(n+1)-b(n+1)-xz_n)/den, K);

        c(n+2) = vpa( ...
            (xy_n-beta*c(n+1))/den, K);
    end

    % Horner法计算 sum_{n=0}^M coefficient_n*h^n
    next = [a(M+1); b(M+1); c(M+1)];

    for n = M-1:-1:0
        next = vpa(h*next + [a(n+1); b(n+1); c(n+1)], K);
    end
end

注意：

y₀、z₀ 精确为 1，表示不对它们施加初始扰动；后续 y(t)、z(t) 仍正常演化，也存在截断与舍入误差。
K = 200 是工作精度，不保证结果有 200 位正确数字。混沌轨道的可信时间需要通过提高 M、提高 K，以及减小步长进行收敛检查。
100 阶、200 位的符号工具箱运算较慢，因此默认 T = 1。微扰极小时，两条轨道在双精度绘图中重合是正常现象，应查看高精度计算后的差距曲线。

lorenze