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

Burgers方程的守恒格式求解

问题描述

我们对下述Burgers方程的初边值问题 { u t + ( u 2 / 2 ) x = 0 − 1 ≤ x ≤ 1 , t > 0 u ( x , 0 ) = u 0 ( x ) = s i n ( π x + π ) p e r i o d i c    B . C .    f o r    x \left \{ \begin{aligned} & u_t + (u^2/2)_x = 0 & -1\leq x \leq 1,t>0\\ & u(x,0) = u_0(x) = sin(\pi x + \pi) & \\ &periodic ~~B.C. ~~for ~~x \end{aligned} \right. ut+(u2/2)x=0u(x,0)=u0(x)=sin(πx+π)periodic  B.C.  for  x1x1,t>0

分别使用不同的守恒格式求解,包括Lax-Friedrichs、Roe(迎风)、Engquist-Osher、Godunov格式和Lax-Wendorff格式。并且要求:

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

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

关键算法格式

精确解

对于如([q])所示的黎曼问题的精确解,我们用特征线方法来求解精确值。过某个初值点 ( x 0 , 0 ) (x_0,0) (x0,0)的特征线是一根直线,它的斜率为 u 0 ( x ) u_0(x) u0(x),那么给定一点 ( x 1 , t 1 ) (x_1,t_1) (x1,t1),我们知道它满足,
x 1 − x 0 t 1 = u 0 ( x 0 ) = s i n ( π x 0 + π ) . \frac{x_1-x_0}{t_1} = u_0(x_0) = sin(\pi x_0 + \pi). t1x1x0=u0(x0)=sin(πx0+π).
求得 x 0 x_0 x0后,我们可以得到在 ( x 1 , t 1 ) (x_1,t_1) (x1,t1)点出的值为 u 0 ( x 0 ) u_0(x_0) u0(x0)
需要注意的是,matlab的solve函数主要是用于多项式求根,对于非线性方程来说,也可以求,但是往往只给出一个最小解。在matlab里面solve命令主要是用来求解代数方程(即多项式)的解,但是也不是说其它方程一个也不能解,不过求解非代数方程的能力相当有限,通常只能给出很特殊的实数解。(该问题给出的方程就是典型的超越方程,非代数方程),在t接近0.5时(比如说t=0.4),x在原点附近(比如x=0.01),就会解错。所以改用fzero函数。

常用守恒格式

常用的守恒格式的写法为: u j n + 1 = u j n − λ ( f ^ j + 1 2 n − f ^ j − 1 2 n ) u_j^{n+1} = u_j^n - \lambda (\hat {f}_{j+\frac{1}{2}}^n-\hat {f}_{j-\frac{1}{2}}^n) ujn+1=ujnλ(f^j+21nf^j21n)
这里的 λ = △ t △ x \lambda = \frac{\triangle t}{\triangle x} λ=xt △ t \triangle t t △ x \triangle x x分别为时间步长和空间步长。
下面提到的几种格式,也是基于这样的一个守恒格式写法,只是数值通量函数 f ^ \hat f f^有所不同。

Lax-Friedrichs格式

Lax-Friedrichs格式的数值通量函数表达式为:
f j + 1 2 n = 1 2 [ f ( u j n ) + f ( u j + 1 n ) − α n ( u j + 1 n − u j n ) ] f_{j+\frac{1}{2}}^n = \frac{1}{2}[f(u_j^n)+f(u_{j+1}^n)-\alpha _n (u_{j+1}^n - u_j^n)] fj+21n=21[f(ujn)+f(uj+1n)αn(uj+1nujn)]
其中,
α n = max ⁡ u n ∣ f ′ ( u ) ∣ \alpha_n = \mathop {\max }\limits_{\rm{u^n}} |f'(u)| αn=unmaxf(u)

Local Lax-Friedrichs格式

Local Lax-Friedrichs格式的数值通量函数表达式为:
f j + 1 2 n = 1 2 [ f ( u j n ) + f ( u j + 1 n ) − α j + 1 2 n ( u j + 1 n − u j n ) ] f_{j+\frac{1}{2}}^n = \frac{1}{2}[f(u_j^n)+f(u_{j+1}^n)-\alpha _{j+\frac{1}{2}}^n (u_{j+1}^n - u_j^n)] fj+21n=21[f(ujn)+f(uj+1n)αj+21n(uj+1nujn)]
其中, α j + 1 2 n = max ⁡ u ∈ ( u j n , u j + 1 n ) ∣ f ′ ( u ) ∣ \alpha _{j+\frac{1}{2}}^n = \mathop {\max }\limits_{u \in (u_j^{n},u_{j+1}^{n})} |f'(u)| αj+21n=u(ujn,uj+1n)maxf(u)

Roe格式

Roe(迎风)格式的数值通量函数为(为了方便省略了第 n n n时间层上标,下如是):
f ^ j + 1 2 = { f ( u j ) f ( u j + 1 ) − f ( u j ) u j + 1 − u j ≥ 0 f ( u j + 1 ) f ( u j + 1 ) − f ( u j ) u j + 1 − u j < 0 \hat f _{j+\frac{1}{2}} = \left \{ \begin{aligned} & f(u_j) & \frac{f(u_{j+1})-f(u_{j})}{u_{j+1}-u_j} \geq 0\\ & f(u_{j+1}) & \frac{f(u_{j+1})-f(u_{j})}{u_{j+1}-u_j} < 0 \end{aligned} \right. f^j+21=f(uj)f(uj+1)uj+1ujf(uj+1)f(uj)0uj+1ujf(uj+1)f(uj)<0

Engquist-Osher格式

数值通量函数的表达式为: f ^ j + 1 2 = f + ( u j ) + f − ( u j + 1 ) \hat f_{j+\frac{1}{2}} = f^+(u_j) + f^-(u_{j+1}) f^j+21=f+(uj)+f(uj+1) 其中,
f + ( u ) = ∫ 0 u max ( f ′ ( v ) , 0 ) d v + f ( 0 ) f − ( u ) = ∫ 0 u min ( f ′ ( v ) , 0 ) d v \begin{aligned} & f^+(u) =& &\int _{0}^{u} \text{max}(f'(v),0)dv + f(0)\\ & f^-(u) =&& \int _{0}^{u} \text{min}(f'(v),0)dv \end{aligned} f+(u)=f(u)=0umax(f(v),0)dv+f(0)0umin(f(v),0)dv 对于Burgers方程而言,有 f + ( u ) = = { f ( u ) u ≥ 0 0 u < 0 f^+(u) = = \left \{ \begin{aligned} & f(u) & u \geq 0\\ & 0 & u < 0 \end{aligned} \right. f+(u)=={f(u)0u0u<0 对于Burgers方程而言,有 f − ( u ) = = { f ( u ) u ≤ 0 0 u > 0 f^-(u) = = \left \{ \begin{aligned} & f(u) & u \leq0\\ & 0 & u > 0 \end{aligned} \right. f(u)=={f(u)0u0u>0

Godunov格式

对于Godunov格式,它的数值通量表达式为: f ^ j + 1 2 = { max ⁡ u j ≤ u ≤ u j + 1 f ( u ) u j ≤ u j + 1 min ⁡ u j ≤ u ≤ u j + 1 f ( u ) u j ≥ u j + 1 \hat f_{j+\frac{1}{2}} = \left \{ \begin{aligned} & \mathop {\max}\limits_{u_j \leq u \leq u_{j+1}}f(u) & u_j \leq u_{j+1}\\ & \mathop {\min}\limits_{u_j \leq u \leq u_{j+1}} f(u) & u_j \geq u_{j+1} \end{aligned} \right. f^j+21=ujuuj+1maxf(u)ujuuj+1minf(u)ujuj+1ujuj+1 对于Burgers方程而言,有 f ^ j + 1 2 = { f ( u j + 1 ) u j ≤ 0 , u j + 1 ≤ 0 f ( u j ) u j ≥ 0 , u j + 1 ≥ 0 f ( max ( ∣ u j ∣ , ∣ u j + 1 ∣ ) ) u j > 0 , u j + 1 < 0 0 u j < 0 , u j + 1 > 0 \hat f_{j+\frac{1}{2}} = \left \{ \begin{aligned} & f(u_{j+1}) & u_j \leq 0 ,u_{j+1} \leq 0\\ & f(u_{j}) & u_j \geq 0 ,u_{j+1} \geq 0\\ & f(\text{max}(|u_{j}|,|u_{j+1}|)) & u_j > 0 ,u_{j+1} < 0\\ & 0 & u_j < 0 ,u_{j+1} > 0 \end{aligned} \right. f^j+21=f(uj+1)f(uj)f(max(uj,uj+1))0uj0,uj+10uj0,uj+10uj>0,uj+1<0uj<0,uj+1>0

Lax-Wendroff格式

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

f ^ j + 1 2 = 1 2 [ f ( u j ) + f ( u j + 1 ) − λ f ′ ( u j + 1 2 ) ( f ( u j + 1 ) − f ( u j ) ) ] \hat f _{j+\frac{1}{2}} = \frac{1}{2}[f(u_j)+f(u_{j+1})-\lambda f'(u_{j+\frac{1}{2}})(f(u_{j+1})-f(u_{j}))] f^j+21=21[f(uj)+f(uj+1)λf(uj+21)(f(uj+1)f(uj))]

这里的 λ \lambda λ含义同上, u j + 1 2 u_{j+\frac{1}{2}} uj+21的值可以通过插值得到,我选的是线性插值,即, u j + 1 2 = 1 2 ( u j + u j + 1 ) u_{j+\frac{1}{2}} = \frac{1}{2}(u_j+u_{j+1}) uj+21=21(uj+uj+1)

数值实验与结果

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

Lax-Friedrichs格式

这里写图片描述
Local Lax-Friedrichs格式

这里写图片描述
Roe格式

这里写图片描述
这里写图片描述
Engquist-Osher格式

这里写图片描述

Godunov格式

这里写图片描述
这里写图片描述
Lax-Wendroff格式

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

这里写图片描述

程序源码

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

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

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

这里写图片描述

已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 点我我会动 设计师:上身试试 返回首页