双曲守恒律的ENO格式和WENO格式.zip

双曲守恒律的ENO格式和WENO格式的MATLAB代码，求解的是Burgers方程的初值问题，分别使用有限体积法4阶ENO格式，有限体积法3阶和5阶WENO方法求解。时间方向用三阶TVD Runge-Kutta方法。研究格式在解光滑情况下和有间断的情况下的数值精度并作图。说明文档参考博客https://blog.csdn.net/lusongno1/article/details/79070049。

function res = ENO4_L（wfluxdfluxSdx）
%% Lax-Friedrichs 通量分裂，u是左值，v是右值
a=max（abs（dflux（w）））; v0=0.5*（flux（w）+a*w）; u0=circshift（0.5*（flux（w）-a*w）[0-1]）;

%% Right Flux
hn = v_reconstructe（v0）;
hp = fliplr（v_reconstructe（fliplr（u0）））;
res = （hp-circshift（hp[01]）+hn-circshift（hn[01]））/dx - S（w）;
end
function hn = v_reconstructe（v0）
v = [v0（end-2:end） v0 v0（1:3）];%前后扩充
n = length（v）;k = 4;m = k-1;
DQ_table = zeros（n4）;
DQ_table（:1） = v;
for i=2:4
    current_col = DQ_table（:i-1）;
    DQ_table（1:end-i+1i） = current_col（2:end-i+2） - current_col（1:end-1-i+2）;
end
IS = 1:n;
abs_DQ = abs（DQ_table）;
for i = 4:n-3
for j = 2:m+1
    if（abs_DQ（IS（i）-1j）        IS（i） = IS（i） - 1;
    end
end
end
Crj = [1/4 13/12 -5/12 1/12;
    -1/12 7/12 7/12 -1/12;
    1/12 -5/12 13/12 1/4;
    -1/4 13/12 -23/12 25/12;
    ];

hn = zeros（）;
for i = 4:n-3
    r = i - IS（i）;
    row_ind = r + 1;
    hn（i-3） = dot（Crj（row_ind:）v（IS（i）:IS（i）+k-1））;
end
end





% vmmm = circshift（v[0 3]）;
% vmm = circshift（v[0 2]）;
% vm  = circshift（v[0 1]）;
% vp  = circshift（v[0 -1]）;
% vpp = circshift（v[0 -2]）;
% vppp = circshift（v[0 -3]）;
% 
% IS = 1:len;%模板起始下标
% for i = 1:len
%     
% end
% 
% 
% % 多项式重构右值
% p0n = （2*vmm - 7*vm + 11*v）/6;%i-2+0为起始的模板
% p1n = （ -vm  + 5*v  + 2*vp）/6;%i-2+1为起始的模板
% p2n = （2*v   + 5*vp - vpp ）/6;%i-2+2为起始的模板
% 
% % Smooth Indicators （Beta factors）
% %syms vmm vm v vp vpp
% B0n = 13/12*（vmm-2*vm+v  ）.^2 + 1/4*（vmm-4*vm+3*v）.^2; 
% B1n = 13/12*（vm -2*v +vp ）.^2 + 1/4*（vm-vp）.^2;
% B2n = 13/12*（v  -2*vp+vpp）.^2 + 1/4*（3*v-4*vp+vpp）.^2;
% 
% % Constants
% d0n = 1/10; d1n = 6/10; d2n = 3/10; epsilon = 1e-6;
% 
% % Alpha weights 
% alpha0n = d0n./（epsilon + B0n）.^2;
% alpha1n = d1n./（epsilon + B1n）.^2;
% alpha2n = d2n./（epsilon + B2n）.^2;
% alphasumn = alpha0n + alpha1n + alpha2n;
% 
% % ENO stencils weigths
% w0n = alpha0n./alphasumn;
% w1n = alpha1n./alphasumn;
% w2n = alpha2n./alphasumn;
% 
% % Numerical Flux at cell boundary $u_{i+1/2}^{-}$;
% hn = w0n.*p0n + w1n.*p1n + w2n.*p2n;

%% Left Flux 
% Choose the negative fluxes ‘u‘ to compute the left cell boundary flux:
% $u_{i-1/2}^{+}$ 

% umm = circshift（u[0 2]）;%m为-，p为+
% um  = circshift（u[0 1]）;
% up  = circshift（u[0 -1]）;
% upp = circshift（u[0 -2]）;
% 
% % Polynomials
% p0p = （ -umm + 5*um + 2*u  ）/6;
% p1p = （ 2*um + 5*u  - up   ）/6;
% p2p = （11*u  - 7*up + 2*upp）/6;
% 
% % Smooth Indicators （Beta factors）
% B0p = 13/12*（umm-2*um+u  ）.^2 + 1/4*（umm-4*um+3*u）.^2; 
% B1p = 13/12*（um -2*u +up ）.^2 + 1/4*（um-up）.^2;
% B2p = 13/12*（u  -2*up+upp）.^2 + 1/4*（3*u -4*up+upp）.^2;
% 
% % Constants
% d0p = 3/10; d1p = 6/10; d2p = 1/10; epsilon = 1e-6;
% 
% % Alpha weights 
% alpha0p = d0p./（epsilon + B0p）.^2;
% alpha1p = d1p./（epsilon + B1p）.^2;
% alpha2p = d2p./（epsilon + B2p）.^2;
% alphasump = alpha0p + alpha1p + alpha2p;
% 
% % ENO stencils weigths
% w0p = al

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        3275  2018-01-15 17:50  ENO4_L.m
     文件         590  2018-01-15 10:31  ExSolu.m
     文件         217  2018-01-15 22:11  ExSolu_Refine.m
     文件       15218  2018-01-15 21:56  ExSoluDatau.mat
     文件        1829  2018-01-15 16:26  gif动图绘制.m
     文件        3149  2018-01-15 22:27  Proj3main.m
     文件        1231  2018-01-15 14:14  WENO3_L.m
     文件        1936  2018-01-15 14:46  WENO5_L.m