常用守恒格式求解双曲问题（Lax-Friedrichts、local Lax-F、Roe、Engquist-Osher、Godunov、Lax-Wendroff）

Burgers方程的守恒格式求解（Proj No.1）

陆嵩 计算数学

问题描述

我们对下述Burgers方程的初边值问题 ⎧⎩⎨⎪⎪⎪⎪ut+(u2/2)x=0u(x,0)=u0(x)=sin(πx+π)periodic B.C. for x−1≤x≤1,t>0

分别使用不同的守恒格式求解，包括Lax-Friedrichs、Roe（迎风）、Engquist-Osher、Godunov格式和Lax-Wendorff格式。并且要求：

计算各个格式在t=0.15时的误差和收敛阶。

画出t=0.5时的精确解和数值解的对比图，每种格式单独画一个图。

关键算法格式

精确解

对于如([q])所示的黎曼问题的精确解，我们用特征线方法来求解精确值。过某个初值点(x0,0)的特征线是一根直线，它的斜率为u0(x)，那么给定一点(x1,t1)，我们知道它满足，

x1−x0t1=u0(x0)=sin(πx0+π).

求得x0后，我们可以得到在(x1,t1)点出的值为u0(x0)。

需要注意的是，matlab的solve函数主要是用于多项式求根，对于非线性方程来说，也可以求，但是往往只给出一个最小解。在matlab里面solve命令主要是用来求解代数方程（即多项式）的解,但是也不是说其它方程一个也不能解，不过求解非代数方程的能力相当有限，通常只能给出很特殊的实数解。(该问题给出的方程就是典型的超越方程，非代数方程)，在t接近0.5时（比如说t=0.4），x在原点附近（比如x=0.01），就会解错。所以改用fzero函数。

常用守恒格式

常用的守恒格式的写法为： un+1j=unj−λ(f^nj+12−f^nj−12)

这里的λ=△t△x，△t和△x分别为时间步长和空间步长。

下面提到的几种格式，也是基于这样的一个守恒格式写法，只是数值通量函数f^有所不同。

Lax-Friedrichs格式

Lax-Friedrichs格式的数值通量函数表达式为： fnj+12=12[f(unj)+f(unj+1)−αn(unj+1−unj)]

其中， αn=maxun|f′(u)|

Local Lax-Friedrichs格式

Local Lax-Friedrichs格式的数值通量函数表达式为： fnj+12=12[f(unj)+f(unj+1)−αnj+12(unj+1−unj)]

其中， αnj+12=maxu∈(unj,unj+1)|f′(u)|

Roe格式

Roe（迎风）格式的数值通量函数为（为了方便省略了第n时间层上标，下如是）：

f^j+12=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪f(uj)f(uj+1)f(uj+1)−f(uj)uj+1−uj≥0f(uj+1)−f(uj)uj+1−uj<0

Engquist-Osher格式

数值通量函数的表达式为： f^j+12=f+(uj)+f−(uj+1) 其中，

f+(u)=f−(u)=∫u0max(f′(v),0)dv+f(0)∫u0min(f′(v),0)dv 对于Burgers方程而言，有 f+(u)=={f(u)0u≥0u<0 对于Burgers方程而言，有 f−(u)=={f(u)0u≤0u>0

Godunov格式

对于Godunov格式，它的数值通量表达式为： f^j+12=⎧⎩⎨⎪⎪maxuj≤u≤uj+1f(u)minuj≤u≤uj+1f(u)uj≤uj+1uj≥uj+1 对于Burgers方程而言，有 f^j+12=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪f(uj+1)f(uj)f(max(|uj|,|uj+1|))0uj≤0,uj+1≤0uj≥0,uj+1≥0uj>0,uj+1<0uj<0,uj+1>0

Lax-Wendroff格式

Lax-Wendroff格式的数值通量表达式为：

f^j+12=12[f(uj)+f(uj+1)−λf′(uj+12)(f(uj+1)−f(uj))]

这里的λ含义同上，uj+12的值可以通过插值得到，我选的是线性插值，即，uj+12=12(uj+uj+1)。

数值实验与结果

对以上提到的各个格式，分别做以下工作：1、计算各个格式在t=0.15时的误差和收敛阶。2、画出t=0.5时的精确解和数值解的对比图，每种格式单独画一个图。

Lax-Friedrichs格式

{width=”12cm”}\

Local Lax-Friedrichs格式

{width=”12cm”}\

Roe格式

{width=”12cm”}\

Engquist-Osher格式

{width=”12cm”}\

Godunov格式

{width=”12cm”}\

Lax-Wendroff格式

从表中可以看出，Lax-Wendroff格式格式可能是二阶收敛的。事实上，刚开始时收敛阶可能表现得不那么稳定，但是随着步长的缩短，这个格式收敛阶最终是稳定在二阶的。


{width=”12cm”}\

程序源码

源代码，使用的语言为MAMTLAB2014a（盗版）,Windows 10环境。

% %% 真解绘图
% for t = 0:0.1:0.6
%     x = -1:0.01:1;
%     u = ExSolu(x,t);
%     figure(1);
%     plot(x,u);
%     title(['t=' num2str(t)]);
%     drawnow;
% end
%%

%% 1. Set initial value; run the solver; plot the solution
clc
clear
close all
scheme = 'LaxW';%选择格式 LaxF LocLaxF Roe EqOsher Godunov LaxW
Nt = 5;%设置空间步
t0 = 0; tf = 0.5; dt = (tf - t0)/Nt;
x1 = -1; x2 = 1; r = 0.5; dx = dt/r;%网格比设置为0.5
xx = x1:dx:x2;%空间划分
%u0 = ExSolu(xx, 0);%求解初值，离散初值满足真解,也可以取点生成
u0 = -sin(pi*xx);
%plot(xx,u0,'color','r');hold on;ezplot('-sin(pi*x)',xx)
uex = ExSolu(xx, tf);%tf时间后的真值
u = eval([scheme '(u0, xx, dx, dt, Nt)']);
plot(xx, uex, 'r-', xx, u, 'g-');
ymin = min(u) - 0.1; ymax = max(u) + 0.1;
axis([x1, x2, ymin, ymax]);
xlabel('x'); ylabel('u');
legend('Exact Solution', scheme, 'Location', 'SouthWest');
title([scheme '  scheme:  t=' num2str(tf)]);

%% 误差计算
ne = 5;%次数
E = zeros(ne,1); Order=zeros(ne,1);%误差和收敛阶初始化
for in = 1:ne
Nt = 30*2^(in-1);%时间步长两倍增长
t0 = 0; tf = 0.15; dt = (tf - t0)/Nt;
x1 = -1; x2 = 1; r = 0.5; dx = dt/r;
xx = x1:dx:x2;
%u0 = ExSolu(xx, 0);
u0 = -sin(pi*xx);
uex = ExSolu(xx, tf);
%u = LaxF(u0, xx, dx, dt, Nt);
u = eval([scheme '(u0, xx, dx, dt, Nt)']);
E(in) = sum(abs(u-uex)) * dx;%计算L1误差
if in>1
Order(in) = log(E(in-1)/E(in))/log(2);
end
end

function f = fun( u )
%UNTITLED3 此处显示有关此函数的摘要
%   此处显示详细说明
f = 1/2*u.^2;
end

function u = ExSolu( x, t )
%计算精确解，这里x支持向量输入，输出多个精确解
syms x0;
%u = zeros(size(x));%初始化u
u0 = sin(pi*x0+pi);
f = (x - x0) - t*u0;
% ezplot(f);
% hold on;
% syms y;
% ezplot(0*y);
x0_ex = zeros(1,length(x));
for i = 1:length(x)
% figure(2);
% ezplot(f(i));
% drawnow;
%x0_ex(i) = double(solve(f(i)));
f_fun = matlabFunction(f(i));
if x(i) > 0
x0_ex(i) = fzero(f_fun,1);
else
x0_ex(i) = fzero(f_fun,-1);
end
if(abs(x0_ex)>1)
error('计算范围超出初值');
end
%x0_ex(i) = eval(x0_ex(i));
end
u = double(subs(u0,x0,x0_ex));%真解计算
end

function u0 = LaxF(u0, xx, dx, dt, Nt)
% LaxF(u0, xx, dx, dt Nt): u0-initial values; xx-coordinates;
syms u;lambda = dt/dx;
nx = length(u0);%空间节点数
%f_hat = zeros(1,nx-1);
f_diff = diff (fun(u),u);

for it = 1:Nt
t = it*dt;%时间标记
alpha = max(abs(double(subs(f_diff,u,u0))));%求解alpha，此写法便于通量函数修改后也能使用
f_hat = 1/2*(  fun(u0(2:end))+fun(u0(1:end-1)) - alpha*(u0(2:end)-u0(1:end-1)) );
un = zeros(1,nx);
un(2:end-1) = u0(2:end-1) - lambda*(f_hat(2:end) - f_hat(1:end-1));
un(1) = ExSolu(-1,t);%使用真解来求边界条件，其实是周期边界都为0
un(nx) = ExSolu(1,t);
u0 = un;
end

function u0 = LocLaxF(u0, xx, dx, dt, Nt)
% LaxF(u0, xx, dx, dt Nt): u0-initial values; xx-coordinates;
syms u;lambda = dt/dx;
nx = length(u0);%空间节点数
%f_hat = zeros(1,nx-1);
f_diff = diff (fun(u),u);

for it = 1:Nt
t = it*dt;%时间标记
f_diff_abs = abs(double(subs(f_diff,u,u0)));
alpha = max(f_diff_abs(1:end-1),f_diff_abs(2:end));%求解alpha，此写法便于通量函数修改后也能使用
f_hat = 1/2*(  fun(u0(2:end))+fun(u0(1:end-1)) - alpha.*(u0(2:end)-u0(1:end-1)) );
un = zeros(1,nx);
un(2:end-1) = u0(2:end-1) - lambda*(f_hat(2:end) - f_hat(1:end-1));
un(1) = ExSolu(-1,t);%使用真解来求边界条件，其实是周期边界都为0
un(nx) = ExSolu(1,t);
u0 = un;
end

function u0 = Roe(u0, xx, dx, dt, Nt)
lambda = dt/dx;
nx = length(u0);%空间节点数
for it = 1:Nt
t = it*dt;%时间标记
a = (fun(u0(2:end))-fun(u0(1:end-1)))./(u0(2:end) - u0(1:end-1));
ind_pick = (1:nx-1) + double((a < 0));
f_hat = fun(u0(ind_pick));
un = zeros(1,nx);
un(2:end-1) = u0(2:end-1) - lambda*(f_hat(2:end) - f_hat(1:end-1));
un(1) = ExSolu(-1,t);%使用真解来求边界条件，其实是周期边界都为0
un(nx) = ExSolu(1,t);
u0 = un;
end

function u0 = EqOsher(u0, xx, dx, dt, Nt)
% LaxF(u0, xx, dx, dt Nt): u0-initial values; xx-coordinates;
syms u;lambda = dt/dx;
nx = length(u0);%空间节点数
%f_hat = zeros(1,nx-1);
%f_diff = diff (fun(u),u);
%fdf_abs_int = int(f_diff,u);
%fun_dfint = matlabFunction(fdf_abs_int);

for it = 1:Nt
t = it*dt;%时间标记
f_plus = zeros(1,nx);
f_minus = zeros(1,nx);
index_plus = find(u0 > 0);
index_minus = find(u0 < 0);
f_plus(index_plus) = 1/2*u0(index_plus).^2;
f_minus(index_minus) = 1/2*u0(index_minus).^2;
%这里利用了burgers方程的特殊性，求到了其原函数。一般双曲守恒方程，这个部分得改
%导函数的正负部截断的原函数不一定能求得，需要利用数值积分手段求积分
f_hat = f_plus(1:end-1) + f_minus(2:end);
% f_p = zeros(1,nx);
% f_m = zeros(1,nx);
% for i = 1:nx
%     if u0(i)>0
%         f_p(i) = 1/2*(u0(i))^2;
%     else
%         f_m(i) = 1/2*(u0(i))^2;
%     end
% end
% diff_f_max = max(dfun(u0),0);
% dfmax_ave = (diff_f_max(1:end-1) + diff_f_max(2:end))/2;
% cumsum_dfmax_ave = cumsum(dfmax_ave);
% int_temp = cumsum_dfmax_ave*dx;
% f_plus = [0 int_temp(1:end-1)] + fun(0);
% diff_f_min = min(dfun(u0),0);
% dfmin_ave = (diff_f_min(1:end-1) + diff_f_min(2:end))/2;
% cumsum_dfmin_ave = cumsum(dfmin_ave);
% f_minus = cumsum_dfmin_ave*dx;
% f_hat = f_plus + f_minus;

% int_absdf = fun_dfint(u0(2:end)) - fun_dfint(u0(1:end-1));
% f_hat = 1/2*(fun(u0(2:end))+fun(u0(1:end-1))) - abs(1/2*(int_absdf)) ;
un = zeros(1,nx);
un(2:end-1) = u0(2:end-1) - lambda*(f_hat(2:end) - f_hat(1:end-1));
un(1) = ExSolu(-1,t);%使用真解来求边界条件，其实是周期边界都为0
un(nx) = ExSolu(1,t);
u0 = un;
end

function u0 = Godunov(u0, xx, dx, dt, Nt)
% Godunov(u0, xx, dx, dt Nt): u0-initial values; xx-coordinates;
lambda = dt/dx;
nx = length(u0);%空间节点数
for it = 1:Nt
t = it*dt;%时间标记

f_hat = zeros(1,nx-1);
% flag = double(u0(1:end-1) > u0(2:end));
% sign = (-1).^flag;
% for i = 1:nx-1
% f_hat(i) = sign(i)*fun(fminbnd(@(u)(sign(i)*fun(u)),u0(i+flag(i)),u0(i+1-flag(i))));
% end
% for i = 1:nx-1
%     if(u0(i)<=u0(i+1))
%         f_hat_com(i) = fun(fminbnd(@fun,u0(i),u0(i+1)));
%     else
%         f_hat_com(i) = -fun(fminbnd(@(u)(-1)*fun(u),u0(i+1),u0(i)));
%     end
%
% end
% f_hat_com - f_hat
vpa = 0;
for i = 1:nx-1
if(u0(i)<=-vpa&&u0(i+1)<=-vpa)
f_hat(i) = fun(u0(i+1));
else if(u0(i)>=vpa&&u0(i+1)>=vpa)
f_hat(i) = fun(u0(i));
else if (u0(i)>vpa&&u0(i+1)<-vpa)
f_hat(i) = fun(     max(     abs(u0(i)),abs(u0(i+1))     )     );
end
end
end
end
un = zeros(1,nx);
un(2:end-1) = u0(2:end-1) - lambda*(f_hat(2:end) - f_hat(1:end-1));
un(1) = ExSolu(-1,t);%使用真解来求边界条件，其实是周期边界都为0
un(nx) = ExSolu(1,t);
u0 = un;
end

function u0 = LaxW(u0, xx, dx, dt, Nt)
% LaxF(u0, xx, dx, dt Nt): u0-initial values; xx-coordinates;
syms u;lambda = dt/dx;
nx = length(u0);%空间节点数
%f_hat = zeros(1,nx-1);
f_diff = diff (fun(u),u);
dfun = matlabFunction(f_diff);

for it = 1:Nt
t = it*dt;%时间标记
u0_bar = ( u0(1:end-1)+u0(2:end)) /2.0 ;
f_hat = 1/2* (     fun(u0(2:end)) + fun(u0(1:end-1)) - lambda*dfun(u0_bar).*...
(fun(u0(2:end))-fun(u0(1:end-1)))    );
un = zeros(1,nx);
un(2:end-1) = u0(2:end-1) - lambda*(f_hat(2:end) - f_hat(1:end-1));
un(1) = ExSolu(-1,t);%使用真解来求边界条件，其实是周期边界都为0
un(nx) = ExSolu(1,t);
u0 = un;
end

一点小感悟：

一阶收敛的数值格式是很容易实现的，即使在编程过程中，那个小地方写错了，也可能会达到一阶。也就是说，一阶格式有一阶收敛程序不一定是对的，但是没有一阶收敛，那么肯定是错的。

试图将这些程序写成通用的解决一阶标量双曲守恒方程的形式，也就是说只是通量函数

f=12u2

改变，只要在fun函数脚本中修改通量函数，那么这些格式依然能用。有一个小问题在于，写成通用的，多了一些求导等不不要的计算，加大了运算量。在ENO和Godnov格式中，没有写成通用格式。主要考虑到：1、ENO格式中积分的计算，若对于一般的通量函数f，正负部截断函数的原函数不一定能求，那么这是可能要用到一些数值积分手段来求近似积分，这对Burgers方程是多余的，引入了额外误差。2、Godunov格式中，若考虑Matlab内置的fminbnd函数来求最大值最小值，无疑引入了多余的误差，以及大大降低了求解速度，所以，这里针对Burgers函数的特殊性，针对性地求解了相对应的数值通量函数。不具普适性和泛化能力。

在求解真解时，利用了特征线。在求解非线性方程根时，不能用solve函数，而应该fzero函数。因为solve函数用于求解多项式的零点尚可，当它在求解稍微复杂一点的非线性方程的零点时，若方程有多个零点，它往往只能给你找个一个零点，而这个零点未必是我们说要的那个，这样就容易出问题。

一个高阶收敛的格式，刚开始的时候收敛阶可能表现得不大稳定，比如上面提到Lax-Wendroff格式，它的收敛阶一会1.6一会儿2.3，但是随着网格的加密，它往往会稳定到它真正的收敛阶。如下所示：