间断初值双曲守恒问题的Lax-Friedrichs和后向欧拉数值解法

间断初值双曲守恒问题的Lax-Friedrichs和后向欧拉数值解法

我们用求解PDE常用的Backwand Euler和Lax-Friedrichs数值方法来求解如下双曲问题:

首先回顾一下双曲问题中一些经典的有限差分格式：

- 该黎曼问题的精确解为：

计算精确解Matlab代码：

function u = ExSolu( x, t )%计算精确解，这里x支持向量输入，输出多个精确解a = 1;%设置a为1u = zeros(size(x));%初始化uu((x-a*t)<=0) = 1;%真解计算end

- Lax-Friedrichs数值求解格式和代码如下：

function u = LaxF(u0, xx, dx, dt, Nt)% LaxF(u0, xx, dx, dt Nt): u0-initial values; xx-coordinates;nx = length(u0);%空间节点数r = dt/dx; a = 1; ar = a*r;u = u0; un = zeros(size(u));%要计算的值初始化（新一层）for it =1:Nt% for j=2:nx-1 %内部节点遍历% un(j) = (1-ar)/2 * u(j+1) + (1+ar)/2 * u(j-1);% end%上述循环可以替换为以下，以节省时间jj=2:nx-1;un(jj) = (1-ar)/2 * u(jj+1) + (1+ar)/2 * u(jj-1);t = it*dt;%时间标记un(1) = 1;%左边界条件un(nx) = 0;%右边界条件u = un;%只保留两层的值就够了end

- 隐式欧拉格式：

写成三对角如下：

相应的Matlab代码：

function u = BwdE(u0, xx, dx, dt, Nt)%u0为初值，xx为空间划分，dx为空间步长，dt为时间步长，Nt为时间点数nx = length(u0);%空间节点数r = dt/dx; a = 1; ar = a*r;d = ones(nx-2, 1); d1 = -ar/2 * d; d2 = ar/2 * d;A = spdiags([d1 d d2], -1:1, nx-2, nx-2);%使用spdiags生成（nx-2）*（nx-2）三对角u = u0; ur = u(:);for it = 1:Nt% set boundary values;t = it*dt;%时间标度，就是说计算到哪了u(1) = ExSolu(xx(1), t);%计算精确解的边界条件，这和数值解的边界条件是一致的u(nx) = ExSolu(xx(nx), t);ur(2:nx-1) = u(2:nx-1);%计算右端项ur(2) = ur(2) + ar/2 * u(1);%将左边加的那一项挪到右边来求解ur(nx-1) = ur(nx-1) - ar/2 *u(nx);u(2:nx-1) = A\ur(2:nx-1);endend

- 误差计算

我们可以计算两种格式的误差和收敛阶，我们这里使用L1误差，i.e.

- 主函数：

%% 1. Set initial value; run the solver; plot the solutionclcclearscheme = 'LaxF';%选择后向欧拉格式Nt = 200;%设置空间步t0 = 0; tf = 1; dt = (tf - t0)/Nt;x1 = -2; x2 = 2; r = 0.5; dx = dt/r;%网格比设置为0.5xx = x1:dx:x2;%空间划分u0 = ExSolu(xx, 0);%求解初值，离散初值满足真解uex = ExSolu(xx, tf);%tf时间后的真值%u = LaxF(u0, xx, dx, dt, Nt);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('Lax-Friedrichs scheme: t=1');%% 误差计算ne = 5;%次数E = zeros(ne,1); Order=zeros(ne,1);%误差和收敛阶初始化for in = 1:neNt = 100*2^(in-1);%时间步长两倍增长t0 = 0; tf = 1; dt = (tf - t0)/Nt;x1 = -2; x2 = 2; r = 0.5; dx = dt/r;xx = x1:dx:x2;u0 = ExSolu(xx, 0);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>1Order(in) = log(E(in-1)/E(in))/log(2);endend

通过计算收敛阶，我们可以知道，这两个非守恒的格式收敛阶皆为0.5。

tip：1、在使用matlab内置求解方程组时，用A\b比inv(A)*b速度要快。2、在matlab中，尽量避免用循环，而改用矩阵和向量运算。3、spdiags函数可以生成三对角。4、使用eval函数，将字符串作为命令运行，对于文件操作很有效。5、我们用的Matlab都是盗版的，如果不想用盗版的，可以用Octave替代，它是开源的。