%% =========================================================
%  一维波动方程：二阶中心差分稳定性教学演示（修正版）
%
%  方程：
%      u_tt = c^2 u_xx
%
%  改进：
%  1. 边界条件改为 Neumann: ux = 0 （自由端反射）
%  2. 波到边界后反射但不反号
%  3. 对比 C<1 与 C>1 两种情况
%  4. 首次碰壁前与自由空间解析解对比
%% =========================================================

clear; clc; close all;

%% ---------------- 基本参数 ----------------
L  = 1.0;               % 空间长度
Nx = 201;               % 网格点数
x  = linspace(0, L, Nx);
dx = x(2) - x(1);

c = 1.0;                % 波速
T_end = 1.4;            % 总时间（故意取得长一些，便于看到反射）

x0 = 0.30;              % 初始脉冲中心

%% ---------------- 初始条件 ----------------
% 初始位移
u0 = exp(-400*(x - x0).^2);

% 初始速度
v0 = zeros(size(x));

% 自由空间初始函数（用于首次碰壁前的精确解）
f = @(xx) exp(-400*(xx - x0).^2);

% 首次碰到左右边界的时间
t_hit_left  = x0 / c;
t_hit_right = (L - x0) / c;
t_hit_first = min(t_hit_left, t_hit_right);

%% ---------------- 两组 Courant 数 ----------------
C_list = [0.9, 1.1];

for icase = 1:length(C_list)

    C  = C_list(icase);
    dt = C * dx / c;
    Nt = floor(T_end / dt);
    t  = (0:Nt) * dt;

    fprintf('========================================\n');
    fprintf('Case %d: C = %.2f, dt = %.5f, Nt = %d\n', icase, C, dt, Nt);

    %% =====================================================
    % 初始化
    %% =====================================================
    u_nm1 = u0;                 % n-1层
    u_n   = zeros(size(x));     % n层
    u_np1 = zeros(size(x));     % n+1层

    % -----------------------------------------------------
    % 首步推进：
    % u_i^1 = u_i^0 + dt*v_i^0 + 0.5*C^2*(u_{i+1}^0 - 2u_i^0 + u_{i-1}^0)
    %
    % 对于 Neumann 边界 ux=0:
    % 左端 ghost 点满足 u_0 = u_2
    % 右端 ghost 点满足 u_{N+1} = u_{N-1}
    %
    % 所以边界二阶导数近似：
    % 左端: (u_2 - 2u_1 + u_0)/dx^2 = 2(u_2-u_1)/dx^2
    % 右端: (u_{N-1} - 2u_N + u_{N+1})/dx^2 = 2(u_{N-1}-u_N)/dx^2
    % -----------------------------------------------------
    u_1 = u0;

    % 内点
    for i = 2:Nx-1
        u_1(i) = u0(i) + dt*v0(i) ...
               + 0.5*C^2*(u0(i+1) - 2*u0(i) + u0(i-1));
    end

    % 左右边界（Neumann）
    u_1(1) = u0(1) + dt*v0(1) + C^2*(u0(2) - u0(1));
    u_1(end) = u0(end) + dt*v0(end) + C^2*(u0(end-1) - u0(end));

    u_nm1 = u0;
    u_n   = u_1;

    %% =====================================================
    % 历史量保存（用于 x-t 图）
    %% =====================================================
    Uhist = zeros(Nx, Nt+1);
    Uhist(:,1) = u0(:);
    if Nt >= 1
        Uhist(:,2) = u_1(:);
    end

    amp_hist = zeros(1, Nt+1);
    amp_hist(1) = max(abs(u0));
    if Nt >= 1
        amp_hist(2) = max(abs(u_1));
    end

    %% =====================================================
    % 画图初始化
    %% =====================================================
    fig = figure('Color','w','Position',[80 50 1450 860]);

    % ---------- 子图1：当前波形 ----------
    subplot(2,2,1);
    h_num = plot(x, u0, 'b-', 'LineWidth', 2); hold on;
    h_exact = plot(x, u0, 'r--', 'LineWidth', 1.5);
    grid on;
    xlim([0 L]);
    ylim([-1.2 1.2]);
    xlabel('x');
    ylabel('u(x,t)');
    title(sprintf('当前波形（Neumann边界, C = %.2f）', C));
    legend('数值解','自由空间解析解（首次碰壁前有效）','Location','northeast');

    h_char1 = xline(x0, 'k--', 'LineWidth', 1.2);
    h_char2 = xline(x0, 'k--', 'LineWidth', 1.2);

    txt1 = text(0.02, 0.92, '', 'Units','normalized', ...
        'FontSize', 11, 'BackgroundColor','w');

    % ---------- 子图2：误差 ----------
    subplot(2,2,2);
    h_err = plot(x, zeros(size(x)), 'm-', 'LineWidth', 2);
    grid on;
    xlim([0 L]);
    ylim([-0.2 0.2]);
    xlabel('x');
    ylabel('u_{num}-u_{exact}');
    title('误差（首次碰壁前）');

    % ---------- 子图3：x-t 演化 ----------
    subplot(2,2,3);
    h_img = imagesc(x, t(1:min(2,end)), Uhist(:,1:min(2,end))');
    set(gca, 'YDir', 'normal');
    colorbar;
    xlabel('x');
    ylabel('t');
    title('x-t 演化图');
    hold on;

    % 特征线（从初始脉冲中心出发）
    tt_line = linspace(0, T_end, 400);
    xp = x0 + c*tt_line;
    xm = x0 - c*tt_line;

    maskp = (xp >= 0 & xp <= L);
    maskm = (xm >= 0 & xm <= L);

    plot(xp(maskp), tt_line(maskp), 'w--', 'LineWidth', 1.3);
    plot(xm(maskm), tt_line(maskm), 'w--', 'LineWidth', 1.3);

    % ---------- 子图4：最大幅值 ----------
    subplot(2,2,4);
    if Nt >= 1
        h_amp = plot(t(1:2), amp_hist(1:2), 'k-', 'LineWidth', 2);
    else
        h_amp = plot(t(1), amp_hist(1), 'k-', 'LineWidth', 2);
    end
    grid on;
    xlabel('t');
    ylabel('max|u|');
    title('最大幅值历史');

    sgtitle(sprintf('一维波动方程：二阶中心差分 + 自由端反射边界（C = %.2f）', C), ...
        'FontSize', 16, 'FontWeight', 'bold');

    drawnow;

    %% =====================================================
    % 时间推进
    %% =====================================================
    for n = 2:Nt

        tn = t(n+1);

        % ---------------- 内点更新 ----------------
        for i = 2:Nx-1
            u_np1(i) = 2*u_n(i) - u_nm1(i) ...
                     + C^2*(u_n(i+1) - 2*u_n(i) + u_n(i-1));
        end

        % ---------------- Neumann边界更新 ----------------
        % 左边界：u_x = 0  => ghost: u(0)=u(2)
        u_np1(1) = 2*u_n(1) - u_nm1(1) ...
                 + 2*C^2*(u_n(2) - u_n(1));

        % 右边界：u_x = 0  => ghost: u(N+1)=u(N-1)
        u_np1(end) = 2*u_n(end) - u_nm1(end) ...
                   + 2*C^2*(u_n(end-1) - u_n(end));

        % 保存历史
        Uhist(:,n+1) = u_np1(:);
        amp_hist(n+1) = max(abs(u_np1));

        % =================================================
        % 解析解：仅在首次碰壁前采用自由空间 d'Alembert 解
        % =================================================
        if tn <= t_hit_first
            u_exact = 0.5*(f(x - c*tn) + f(x + c*tn));
            err_now = u_np1 - u_exact;
            show_exact = true;
        else
            u_exact = nan(size(x));
            err_now = nan(size(x));
            show_exact = false;
        end

        %% ---------------- 更新图形 ----------------
        % 当前波形
        subplot(2,2,1);
        set(h_num, 'YData', u_np1);

        if show_exact
            set(h_exact, 'YData', u_exact, 'Visible', 'on');
        else
            set(h_exact, 'Visible', 'off');
        end

        % 当前特征线位置
        delete(h_char1);
        delete(h_char2);

        xp_now = x0 + c*tn;
        xm_now = x0 - c*tn;

        if xp_now >= 0 && xp_now <= L
            h_char1 = xline(xp_now, 'k--', 'LineWidth', 1.2);
        else
            h_char1 = xline(-1, 'k--', 'LineWidth', 1.2); % 图外占位
        end

        if xm_now >= 0 && xm_now <= L
            h_char2 = xline(xm_now, 'k--', 'LineWidth', 1.2);
        else
            h_char2 = xline(-1, 'k--', 'LineWidth', 1.2);
        end

        if C <= 1
            theory_str = '理论上稳定';
        else
            theory_str = '理论上不稳定';
        end

        if show_exact
            maxerr = max(abs(err_now));
            info_str = sprintf('t = %.3f   C = %.2f   max error = %.3e   %s', ...
                tn, C, maxerr, theory_str);
        else
            info_str = sprintf('t = %.3f   C = %.2f   已进入边界反射阶段   %s', ...
                tn, C, theory_str);
        end

        set(txt1, 'String', info_str);

        % 误差图
        subplot(2,2,2);
        if show_exact
            set(h_err, 'YData', err_now, 'Visible', 'on');
            err_lim = max(0.05, 1.2*max(abs(err_now)));
            ylim([-err_lim, err_lim]);
            title('误差（首次碰壁前）');
        else
            set(h_err, 'YData', zeros(size(x)), 'Visible', 'off');
            title('误差（反射后未显示自由空间解析解）');
        end

        % x-t 图
        subplot(2,2,3);
        set(h_img, 'XData', x, 'YData', t(1:n+1), 'CData', Uhist(:,1:n+1)');
        axis([0 L 0 T_end]);

        % 最大幅值历史
        subplot(2,2,4);
        set(h_amp, 'XData', t(1:n+1), 'YData', amp_hist(1:n+1));
        xlim([0 T_end]);

        amp_max = max(amp_hist(1:n+1));
        ylim([0, max(1.2, 1.1*amp_max)]);

        drawnow;

        % 若明显发散则提前停止
        if max(abs(u_np1)) > 1e3
            fprintf('C = %.2f 时数值明显发散，提前停止。\n', C);
            break;
        end

        % 滚动推进
        u_nm1 = u_n;
        u_n   = u_np1;
    end

    disp('按任意键继续下一个案例...');
    pause;
end