计算科学 - 一维无粘性汉堡方程的制造解的问题 - 吾爱随笔录

我对非线性一维 PDE 的最简单示例（无粘性）汉堡方程有疑问：

$u_t + uu_x = 0,~~(1)$

可以重写为一些对流方程

$u_t + f(u)_x = 0$

有助焊剂 $f(u) = (\frac{u}{2})^2$ . 我使用半离散方法来解决具有周期性边界条件的 (1)：我在空间中使用逆风，在时间上使用欧拉向前，总体上得到一阶方案。为了测试收敛性，我使用制造解决方案的方法：

对于一些初始条件 $u(x,t)$ , $u(x,t)$ 是一个解决方案

$u_t + uu_x = r(x,t)$ ,

其中残差 $r(x,t)$ 是结果 $u(x,t)$ 插入 (1)。由于逆风需要正的平流速度，而速度由解决定，u，我选择

$u(x,t) = 5 + \sin(2\pi(x-t)), ~x \in [0,1], ~t \in [0,0.5]$

作为初始解决方案。它认为，

$u_t = -2\pi\cos(2\pi(x-t))$ ,

$u_x = 2\pi\cos(2\pi(x-t))$ ,

所以我的残差是 $r(x,t) = 2\pi\cos(2\pi(x-t))(4 + \sin(2\pi(x-t)))$ . 在计算上风之后，残差作为一些源项处理并添加到更新项。

因为 - 每个构造 - 解决方案应该是 $u(x,t)$ 对于每个 $t \in [0,\infty]$ 和 $x \in [0,1]$ ，我一定遗漏了一些明显的东西。或者有限差分和这个方程是否存在一般问题？我有一个旧的有限体积代码供参考，它工作得很好。

我附上了初始解决方案和解决方案的图 $t=0.25$ 并且 - 由于 matlab 尽可能接近伪代码 - matlab 代码：

clc

format long
N   = 20;   % Number of Points
cfl = 0.5; % 
adv = 2.0; % Linear Advection speed
t_start = 0;
t_end   = 0.25;

% EQ type
% type = "linear_advection";
type    = "burgers";
% fd type
fd      = "upwind";
% fd      = "downwind";
% fd      = "central";
% Initial Conditions
IC      = "default";
% IC      = "resi_test_const";

% switch to plot immediately
plot_immediate = true;

% Initial solution and resdiduals
if type == "burgers"
    if IC == "resi_test_const"
        sol = @(x,t) 5 + sin(2*pi*x);
        r   = @(x,t) sin(2*pi*x)*2*pi.*cos(2*pi*x);
    else
        sol = @(x,t) 5 + sin(2*pi*(x-t));
        r   = @(x,t) 2*pi*cos(2*pi*(x-t)).*(4 + sin(2*pi*(x-t)));
    end
elseif type == "linear_advection"
    if IC == "resi_test_const"
        sol = @(x,t) sin(2*pi*x);
        r   = @(x,t) adv*2*pi*cos(2*pi*x);
    else
        sol = @(x,t) sin(2*pi*(x-t));
        r   = @(x,t) zeros(1,length(x));
    end
end
% Flux
if type == "burgers"
    f = @(u) (u./2).^2;
else
    f = @(u) adv*u;
end

dx  = 1/N;
x   = (0:dx:1-dx)+dx/2;
u   = zeros(1,length(x)+2);
u_t = u;
% Initial Solution
u(2:end-1)  = sol(x,t_start);
% Ghost cells
u(1)        = u(end-1);
u(end)      = u(2);
% Initialize flux
fu = u;

figure(1);
% Plot initial conditions
plot(x,u(2:end-1))
title('Initial Solution', 'Interpreter', 'latex') 
xlabel('x')
ylabel('u')

iter = 1;
t = t_start;
while t<t_end
    % Update dt
    if type == "burgers"
        dt = cfl*0.5*dx/max(abs(u(:)));
    else
        dt = cfl*0.5*dx/abs(adv);
    end
    % Update flux
    fu = f(u);
    if (t+dt)>t_end
        dt = t_end - t;
    end
    if fd == "upwind"
        % Upwinding
        for i=2:length(u)-1
            u_t(i) = (fu(i-1)-fu(i))/dx;
        end
    elseif fd == "downwind"
        for i=2:length(u)-1
            u_t(i) = (fu(i)-fu(i+1))/dx;
        end
    elseif fd == "central"
       for i=2:length(u)-1
           u_t(i) = (fu(i-1)-fu(i+1))/(2*dx);
       end
    end
    % Add source terms
    u_t(2:end-1) = u_t(2:end-1) + r(x,t);
    % Ghost cell update
    u_t(1)        = u_t(end-1);
    u_t(end)      = u_t(2);
    % Update u (euler forward)
    u = u + dt*u_t;
    % Update current time and iteration counter
    iter = iter + 1;
    t = t + dt;
    if plot_immediate
        % Draw plot immediately
        figure(2);
        drawnow
        plot(x,u(2:end-1))
        title(['$t = $', num2str(t), ', $n_{\textrm{cells}} = $', ...
            num2str(N)], 'Interpreter', 'latex') 
        xlabel('x')
        ylabel('u')
    end
end

if ~plot_immediate
    figure(2);
    plot(x,u(2:end-1))
    title(['$t = $', num2str(t), ', $n_{\textrm{cells}} = $', ...
        num2str(N)], 'Interpreter', 'latex') 
    xlabel('x')
    ylabel('u')
end

该代码有一些开关可以解决，例如带有残差的线性平流，因此解决方案是 $u(x,t) = \sin(2\pi x)$ . 这工作得很好，所以我很确定我对汉堡方程、残差或制造解决方案的想法有疑问......非常感谢提示！