%% =========================================================
%  Lax-Wendroff 方法教学演示
%  一维线性对流方程：
%
%       u_t + a u_x = 0
%
%  解析解：
%
%       u(x,t) = u0(x - a t)
%
%  采用周期边界条件
%
%  课堂展示：
%  1. Lax-Wendroff 格式
%  2. CFL 数对稳定性的影响
%  3. 数值解与解析解对比
%  4. 数值色散现象
%% =========================================================

clear; clc; close all;

%% ---------------------------------------------------------
%  1. 参数设置
%% ---------------------------------------------------------

L  = 1.0;          % 计算区域长度
Nx = 200;          % 网格点数
a  = 1.0;          % 对流速度

x  = linspace(0, L, Nx+1);
x(end) = [];       % 去掉最后一个点，避免周期点重复

dx = L / Nx;       % 空间步长

CFL = 0.8;         % CFL 数，Lax-Wendroff 稳定条件为 |CFL| <= 1
dt  = CFL * dx / abs(a);

T_end = 1.0;       % 终止时间
Nt = round(T_end / dt);
dt = T_end / Nt;   % 调整 dt，使得刚好到达 T_end
CFL = a * dt / dx; % 实际 CFL 数

fprintf('dx  = %.5f\n', dx);
fprintf('dt  = %.5f\n', dt);
fprintf('CFL = %.5f\n', CFL);

%% ---------------------------------------------------------
%  2. 初始条件
%% ---------------------------------------------------------

% 高斯脉冲
u0 = exp(-100 * (x - 0.3).^2);

% 数值解初始化
u = u0;

%% ---------------------------------------------------------
%  3. Lax-Wendroff 格式
%
%  对方程：
%
%       u_t + a u_x = 0
%
%  Lax-Wendroff 格式为：
%
%       u_i^{n+1}
%       = u_i^n
%       - CFL/2 * (u_{i+1}^n - u_{i-1}^n)
%       + CFL^2/2 * (u_{i+1}^n - 2u_i^n + u_{i-1}^n)
%
%  其中：
%
%       CFL = a dt / dx
%
%% ---------------------------------------------------------

figure('Color','w');

for n = 1:Nt

    % -----------------------------------------------------
    % 周期边界下的左右邻点
    %
    % ip 表示 i+1
    % im 表示 i-1
    %
    % circshift(u,-1)：整体左移，相当于 u_{i+1}
    % circshift(u, 1)：整体右移，相当于 u_{i-1}
    % -----------------------------------------------------

    u_ip = circshift(u, -1);
    u_im = circshift(u,  1);

    % -----------------------------------------------------
    % Lax-Wendroff 更新
    % -----------------------------------------------------

    u_new = u ...
          - 0.5 * CFL   * (u_ip - u_im) ...
          + 0.5 * CFL^2 * (u_ip - 2*u + u_im);

    u = u_new;

    %% -----------------------------------------------------
    %  每隔若干步画一次图
    %% -----------------------------------------------------

    if mod(n, 5) == 0 || n == 1 || n == Nt

        t = n * dt;

        % -------------------------------------------------
        % 解析解
        %
        % 由于是周期区域，需要把 x - a t 映射回 [0,L)
        % -------------------------------------------------

        x_shift = mod(x - a*t, L);
        u_exact = exp(-100 * (x_shift - 0.3).^2);

        % 误差
        err = u - u_exact;
        L2err = sqrt(sum(err.^2) / Nx);

        % -------------------------------------------------
        % 绘图
        % -------------------------------------------------

        clf;

        subplot(2,1,1);
        plot(x, u_exact, 'k-', 'LineWidth', 2); hold on;
        plot(x, u, 'ro', 'MarkerSize', 4, 'LineWidth', 1);
        grid on;

        xlabel('x');
        ylabel('u');

        title(sprintf('Lax-Wendroff 方法：t = %.3f, CFL = %.2f, L_2 error = %.3e', ...
              t, CFL, L2err));

        legend('解析解', '数值解', 'Location', 'best');

        ylim([-0.2, 1.2]);

        subplot(2,1,2);
        plot(x, err, 'b-', 'LineWidth', 1.5);
        grid on;

        xlabel('x');
        ylabel('误差');

        title('数值误差 u_{num} - u_{exact}');

        ylim([-0.2, 0.2]);

        drawnow;
    end
end

%% ---------------------------------------------------------
%  4. 最终误差输出
%% ---------------------------------------------------------

t = T_end;
x_shift = mod(x - a*t, L);
u_exact = exp(-100 * (x_shift - 0.3).^2);

err = u - u_exact;

L1err   = sum(abs(err)) / Nx;
L2err   = sqrt(sum(err.^2) / Nx);
Linferr = max(abs(err));

fprintf('\n最终误差：\n');
fprintf('L1   error = %.6e\n', L1err);
fprintf('L2   error = %.6e\n', L2err);
fprintf('Linf error = %.6e\n', Linferr);