带源项对流方程解析解与有限差分解的差异问询

我首次使用Stack Overflow，若内容显示异常或表述不清，请随时告知。我正在求解一个已知解析解的带源项对流方程，采用有限差分法离散求解。结果虽相近但存在差异，求解区域分为四段：第一段解析解与数值解完全重合；第二段吻合度较好，但已出现误差；第三、四段误差情况持续恶化。

解析解计算函数 #

function [C1, C2, C3, C4] = calcul_sol_analy(mesures,params)

    
    % Solution portion 1
    x_p1 = linspace(0, params.xcm, 100);        
    [X1, T1] = meshgrid(x_p1, params.t_values);
    C1 = params.c0(T1 - X1/params.v) .* exp(-params.k * X1 / params.v);   

    % Solution portion 2
    Cm =  interp1(params.date_jours_ch, mesures, params.t_values, 'pchip')';  % transposée 
    x_p2 = linspace(params.xcm, params.xcp, 100);  
    [X2, T2] = meshgrid(x_p2, params.t_values);
    c1 = @(t) (C1(:, end) .* params.debit_s + Cm .* params.debit_m) ./ (params.debit_SM);
    C2 = c1(T2 - (X2-params.xcm)/params.v) .* exp(-params.k * (X2-params.xcm)/ params.v);   
 
    % Solution portion 3
    
    Cp = interp1(params.date_jours_cl, params.mesures_crassphage_ch, params.t_values, 'pchip')'; % transposée 
    x_p3 = linspace(params.xcp, params.xco, 100); 
    [X3, T3] = meshgrid(x_p3, params.t_values);
    c2 = @(t) (C2(:, end) .* params.debit_SM + Cp .* params.debit_cl) ./ (params.debit_SME);
    C3 = c2(T3 - (X3-params.xcp)/params.v) .* exp(-params.k * (X3-params.xcp)/ params.v); 

    % Solution portion 4
   
    Co = interp1(params.date_jours_tr, params.mesures_crassphage_ch, params.t_values, 'pchip')'; % transposée
    x_p4 = linspace(params.xco, params.xct, 100); 
    [X4, T4] = meshgrid(x_p4, params.t_values);
    c3 = @(t) (C3(:, end) .* params.debit_SME + Co .* params.debit_o) ./ (params.debit_SMO);
    C4 = c3(T4 - (X4-params.xco)/params.v) .* exp(-params.k * (X4-params.xco)/ params.v);

end

数值解计算函数 #

function [C1, C2, C3, C4] = calcul_sol_num(mesures, params)

    % Solution portion 1
    C1 = zeros(params.Nt, params.Nx);
    C1(:, 1) = params.c0(params.t_values);
    x_p1 = linspace(0, params.xcm, params.Nx);  
    [X1, T1] = meshgrid(x_p1, params.t_values);
    C_exact = params.c0(T1 - X1 / params.v) .* exp(-params.k * X1 / params.v);
    C1(1, :) = C_exact(1, :);
    
    for n = 1:params.Nt-1
        for i = 2:params.Nx-1
            C1(n + 1, i) = C1(n, i) * (1 - params.k * params.dt - params.v * params.dt / params.dx) + ...
                           C1(n, i-1) * (params.v * params.dt / params.dx);
        end
        C1(n + 1, params.Nx) = C1(n + 1, params.Nx-1);
    end

    % Solution portion 2
    Cm =  interp1(params.date_jours_ch, mesures, params.t_values, 'pchip')';  % transposée 
    x_p2 = linspace(params.xcm, params.xcp, params.Nx);  
    [X2, T2] = meshgrid(x_p2, params.t_values);
    c1 = @(t) (C1(:, end) .* params.debit_s + Cm .* params.debit_m) ./ params.debit_SM;
    C2 = zeros(params.Nt, params.Nx);
    C2(:, 1) = c1(params.t_values);
    C_exact = c1(T2 - (X2-params.xcm) / params.v) .* exp(-params.k * (X2-params.xcm) / params.v);
    C2(1, :) = C_exact(1, :);
    
    for n = 1:params.Nt-1
        for i = 2:params.Nx-1
            C2(n + 1, i) = C2(n, i) * (1 - params.k * params.dt - params.v * params.dt / params.dx) + ...
                           C2(n, i-1) * (params.v * params.dt / params.dx);
        end
        C2(n + 1, params.Nx) = C2(n + 1, params.Nx-1);
    end 

    % Solution portion 3
    Cp =  interp1(params.date_jours_cl, params.mesures_crassphage_ch, params.t_values, 'pchip')';  % transposée 
    x_p3 = linspace(params.xcp, params.xco, params.Nx);  
    [X3, T3] = meshgrid(x_p3, params.t_values);
    c2 = @(t) (C2(:, end) .* params.debit_SM + Cp .* params.debit_cl) ./ params.debit_SME;
    C3 = zeros(params.Nt, params.Nx);
    C3(:, 1) = c2(params.t_values);
    C_exact = c2(T3 - (X3-params.xcp) / params.v) .* exp(-params.k * (X3-params.xcp) / params.v);
    C3(1, :) = C_exact(1, :);
    
    for n = 1:params.Nt-1
        for i = 2:params.Nx-1
            C3(n + 1, i) = C3(n, i) * (1 - params.k * params.dt - params.v * params.dt / params.dx) + ...
                           C3(n, i-1) * (params.v * params.dt / params.dx);
        end
        C3(n + 1, params.Nx) = C3(n + 1, params.Nx-1);
    end 

    % Solution portion 4
    Co = interp1(params.date_jours_tr, params.mesures_crassphage_oi, params.t_values, 'pchip')'; 
    x_p4 = linspace(params.xco, params.xct, params.Nx);  
    [X4, T4] = meshgrid(x_p4, params.t_values);
    c3 = @(t) (C3(:, end) .* params.debit_SME + Co .* params.debit_o) ./ params.debit_SMO;
    C4 = zeros(params.Nt, params.Nx);
    C4(:, 1) = c3(params.t_values);
    C_exact = c3(T4 - (X4-params.xco) / params.v) .* exp(-params.k * (X4-params.xco)  / params.v);
    C4(1, :) = C_exact(1, :);
    
    for n = 1:params.Nt-1
        for i = 2:params.Nx-1
            C4(n + 1, i) = C4(n, i) * (1 - params.k * params.dt - params.v * params.dt / params.dx) + ...
                           C4(n, i-1) * (params.v * params.dt / params.dx);
        end
        C4(n + 1, params.Nx) = C4(n + 1, params.Nx-1);
    end 

end

离散参数 #

Nx = 101;  

Tf = max(date_jours_ch);

dx = L / (Nx-1) ;
dt = dx / v;
x_values = 0:dx:L;              
t_values = 2:dt:Tf;  
Nt=length(t_values);

各段解对比 #

• 第一段解：解析解与数值解完全重合
• 第二段解：吻合度较好，但已出现误差
• 第三段解：误差情况持续恶化
• 第四段解：误差进一步加剧

恳请各位帮忙分析误差逐段恶化的原因，以及如何提升数值解的精度。

内容的提问来源于stack exchange，提问作者thr