% u_t + f(u)_x =0 u(x,0) = 0.5+sin(pi*x)
% v_t + g(u,v)_x =0 v(x,0) = pi*cos(pi*x)
% f(u) = u ,g(u,v) = v Cycle boundary conditions tend=1/(2*pi)
% Exact solution：u(x,t) = u0(x-t)=0.5+sin(pi*(x-t))
clear; clc; close all; format long;%Main program
er = zeros(1,7);
for i = 1:7
n = 10*2^(i-1);%n and so on
er(i) = error(n);
end
for i = 1:6
order = (log(er(i))-log(er(i+1)))/log(2)
end
function f = flux(u)%An expression of the flux function f(u).
f = u;
end
function g = glux(u,v)%An expression of the flux function g(u,v).
g = v;
end
function y = fun(x)%The expression of the initial value of u
y = 0.5 + sin(pi*x);
end
function y = dfun(x)%An expression of the initial value of v
y = pi * cos(pi*x);
end
function y = gun(x)%Exact solution expression on T
y = 0.5 + sin(pi*(x-1/(2*pi)));
end
function [u,v] = hweno5_rk3(n) %Calculate the you and v (mean) of the last layer
global te xa xb dx
xa = 0;
xb = 2;
dx = (xb - xa)/n;%Karculat has Andvi (Gate) Bentrastrell
x = xa:dx:xb;%Half a grid, n+1
cflnum = 0.2;
te = 1/(2*pi);
t = 0;u = zeros(1,n);v = zeros(1,n);
for i = 1:n
u(i) = comsimpson(@fun,x(i),x(i+1),5000) * (1/dx);%The mean of the u-grid for all the whole points of the first layer
v(i) = comsimpson(@dfun,x(i),x(i+1),5000) * (1/dx);%The mean of the v-grid for all the whole points of the first layer
end
%Time-level advancement
while t < te
alpha = 1;
dt = (cflnum*(dx^(5/3)))/alpha;
if t + dt <= te
t = t + dt;
else
dt = te- t;
t = te;
end
[u,v] = rk(u,v,n,dt);
end
end
function [m,n] = HWENO5(a,b,c,d,e,f,dx)%Positive and negative universal, the first three are u_i-1, u_i, u_i+1; The last three are v_i-1, v_i, v_i+1
epsilon = 1e-8;
%m reconstruction:u
p0u = (-7*a + 13*b - 4*dx*d)/6;
p1u = (b + 5*c - 2*dx*f)/6;
p2u = (-a + 5*b + 2*c)/6;
beta0u = (-2*a + 2*b - dx*d)^2 + 13/3*(-a + b - dx*d)^2;
beta1u = (-2*b + 2*c - dx*f)^2 + 13/3*(-b + c - dx*f)^2;
beta2u = 1/4*(-a + c)^2 + 13/12*(-a + 2*b - c)^2;
w0u = 9/(80*(epsilon + beta0u)^2);
w1u = 21/(40*(epsilon + beta1u)^2);
w2u = 29/(80*(epsilon + beta2u)^2);
Su = w0u + w1u + w2u;
w0u = w0u/Su;
w1u = w1u/Su;
w2u = w2u/Su;
m = w0u*p0u + w1u*p1u + w2u*p2u;
%n reconstruction:v
p0v = (8*a - 8*b + 3*dx*d + 7*dx*e)/(2*dx);
p1v = (-4*b + 4*c - dx*e - dx*f)/(2*dx);
p2v = (a - 4*b + 3*c + 2*dx*e)/(4*dx);
beta0v = 4*(3*a - 3*b + dx*d + 2*dx*e)^2 + 39*(2*a - 2*b + dx*d + dx*e)^2;
beta1v = 4*(-3*b + 3*c - 2*dx*e - dx*f)^2 + 39*(2*b - 2*c + dx*e + dx*f)^2;
beta2v = (a -2*b + c)^2 + 39/4*(-a + c - 2*dx*e)^2;
w0v = 1/(18*(epsilon + beta0v)^2);
w1v = 5/(6*(epsilon + beta1v)^2);
w2v = 1/(9*(epsilon + beta2v)^2);
Sv = w0v + w1v + w2v;
w0v = w0v/Sv;
w1v = w1v/Sv;
w2v = w2v/Sv;
n = w0v*p0v + w1v*p1v + w2v*p2v;
end
function [du,dv] = Lu(u,v,n)%Find the right term
global dx
ub = [[0,0],u,[0,0]];%The mean of u after the expansion of the periodic condition
vb = [[0,0],u,[0,0]];%The mean of v after the expansion of the periodic condition
un = zeros(1,n + 1);%The half-point reconstruction value of u-
up = zeros(1,n + 1);%The half-point reconstruction value of u+
vn = zeros(1,n + 1);%The half-point reconstruction value of v-
vp = zeros(1,n + 1);%The half-point reconstruction value of v+
fhat = zeros(1,n + 1);
ghat = zeros(1,n + 1);
du = zeros(1,n);
dv = zeros(1,n);
% Cycle condition expansion
ub(1) = u(n-1); vb(1) = v(n-1);
ub(2) = u(n); vb(2) = v(n);
ub(end-1) = u(1); vb(end-1) = v(1);
ub(end) = u(2); vb(end) = v(2);
for i = 1:n + 1
[un(i),vn(i)] = HWENO5(ub(i),ub(i + 1),ub(i + 2),vb(i),vb(i + 1),vb(i + 2),dx);%-
[up(i),vp(i)] = HWENO5(ub(i + 3),ub(i + 2),ub(i + 1),vb(i + 3),vb(i + 2),vb(i + 1),dx);%+
end
for i = 1:n + 1
alpha = 1;
fhat(i) = 1/2*(flux(un(i)) + flux(up(i)) - alpha *(up(i)-un(i)));%Calculate the LAX flux at all half points J+1/2, and the space layer advances
ghat(i) = 1/2*(glux(un(i),vn(i)) + glux(up(i),vp(i)) - alpha *(vp(i)-vn(i)));
end
for i = 1:n
du(i) = -(1/dx)*(fhat(i + 1) - fhat(i));
dv(i) = -(1/dx)*(ghat(i + 1) - ghat(i));
end
end
function [ux,vx] = rk(u,v,n,dt)%RK-3
[du,dv] = Lu(u,v,n);
u1 = u + (dt)*du; v1 = v + (dt)*dv;
[du,dv] = Lu(u1,v1,n);
u2 =(3/4)*u + (1/4)*u1 + (1/4)*dt*du; v2 =(3/4)*v + (1/4)*v1 + (1/4)*dt*dv;
[du,dv] = Lu(u2,v2,n);
ux = (1/3)*u + (4/6)*u2 + (4/6)*dt*du; vx = (1/3)*v + (4/6)*v2 + (4/6)*dt*dv;
end
function yj = comsimpson(fun,a,b,d)
%The complex Simpson formula is integral, fun is the integrable function, ab is the upper limit under the integral, d is the number of intervals, and the returned y is a numerical result
z1 = feval(fun,a)+feval(fun,b);m = d/2;
h = (b-a)/(2*m);x=a;
z2 = 0;z3 = 0;x2 = 0;x3 = 0;
for k = 2:2:2*m
x2 = x + k*h;z2 = z2 + 2*feval(fun,x2);
end
for k = 3:2:2*m
x3 = x + k*h;z3 = z3 + 4*feval(fun,x3);
end
yj = (z1 + z2 + z3) * h/3;
end
function e = error(n)
global xa xb dx
dx = (xb - xa)/n;
x = xa:dx:xb;
uu = zeros(1,n);%Exact solution of mesh integrals
for i = 1:n
uu(i) = comsimpson(@gun,x(i),x(i+1),5000) * (1/dx);%The mean of the mesh for all the whole points of the first layer
end
[u,v] = hweno5_rk3(n);%u is the numerical solution;
e = sum(abs(u-uu))/(n+1)%Calculation error
end
