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

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

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

这里写图片描述

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

这里写图片描述

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

这里写图片描述

计算精确解Matlab代码:

function u = ExSolu( x, t )
%计算精确解,这里x支持向量输入,输出多个精确解
a = 1;%设置a为1
u = zeros(size(x));%初始化u
u((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);
end
end

- 误差计算

我们可以计算两种格式的误差和收敛阶,我们这里使用 L1 误差,i.e.
这里写图片描述

- 主函数:

%% 1. Set initial value; run the solver; plot the solution
clc
clear
scheme = 'LaxF';%选择后向欧拉格式
Nt = 200;%设置空间步
t0 = 0; tf = 1; dt = (tf - t0)/Nt;
x1 = -2; x2 = 2; r = 0.5; dx = dt/r;%网格比设置为0.5
xx = 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:ne
Nt = 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>1
Order(in) = log(E(in-1)/E(in))/log(2);
end
end

这里写图片描述

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

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

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