一维burgers方程的RKDG有限元解法 matlab代码

以下是一维Burgers方程的RKDG有限元解法的Matlab代码：

clear all;
clc;
L = 2*pi; % 定义空间区间[0,L]
N = 100; % 定义空间网格数
dx = L/N; % 计算空间步长
x = 0:dx:L; % 定义空间网格
nu = 0.1; % 定义粘性系数
dt = 0.001; % 定义时间步长
T = 2*pi; % 定义时间区间[0,T]
Nt = round(T/dt); % 计算时间步数
cfl = 0.8; % 定义CFL数
u = sin(x); % 定义初始条件
u1 = u; % 保存初始条件
Nq = 3; % 定义积分精度
Ng = Nq+1; % 定义高斯积分点数
Np = 2; % 定义多项式次数
K = N-1; % 定义单元数
v = zeros(K+1,Np+1); % 定义测试函数
B = zeros(K+1,Np+1); % 定义基函数
for k=1:K+1 % 循环遍历每个单元
    xl = x(k); % 定义单元左端点
    xr = x(k+1); % 定义单元右端点
    xq = zeros(1,Ng); % 定义高斯积分点
    wq = zeros(1,Ng); % 定义高斯积分权重
    [xq,wq] = gauss(Nq); % 计算高斯积分点和权重
    for i=1:Np+1 % 循环遍历每个基函数
        xi = (xl + xr)/2 + (xl - xr)/2 * xq(i); % 计算基函数中心位置
        B(k,i) = lagrange_basis(xi,xl,xr,i-1); % 计算基函数值
        v(k,i) = B(k,i); % 测试函数取基函数值
    end
end
A = zeros(K+1,Np+1,K+1,Np+1); % 定义刚度矩阵
M = zeros(K+1,Np+1,K+1,Np+1); % 定义质量矩阵
for k1=1:K+1 % 循环遍历每个单元
    for i1=1:Np+1 % 循环遍历每个基函数
        for k2=1:K+1 % 循环遍历每个单元
            for i2=1:Np+1 % 循环遍历每个基函数
                A(k1,i1,k2,i2) = A(k1,i1,k2,i2) + nu/dx * dot(d_basis(x(k1),i1-1),d_basis(x(k2),i2-1)) * dx/2; % 计算刚度矩阵
                M(k1,i1,k2,i2) = M(k1,i1,k2,i2) + dot(B(k1,i1),B(k2,i2)) * dx/2; % 计算质量矩阵
            end
        end
    end
end
A = reshape(reshape(A,(K+1)*(Np+1),(K+1)*(Np+1)),(K+1)*(Np+1),(K+1)*(Np+1)); % 将刚度矩阵转换为向量
M = reshape(reshape(M,(K+1)*(Np+1),(K+1)*(Np+1)),(K+1)*(Np+1),(K+1)*(Np+1)); % 将质量矩阵转换为向量
B = reshape(B,(K+1)*(Np+1),1); % 将基函数转换为向量
v = reshape(v,(K+1)*(Np+1),1); % 将测试函数转换为向量
for n=1:Nt % 循环遍历每个时间步长
    un = u; % 保存当前时间步长下的速度场
    for k=1:K+1 % 循环遍历每个单元
        xl = x(k); % 定义单元左端点
        xr = x(k+1); % 定义单元右端点
        uq = zeros(1,Ng); % 定义速度场插值点
        for i=1:Ng % 循环遍历每个高斯积分点
            xi = (xl + xr)/2 + (xl - xr)/2 * xq(i); % 计算高斯积分点位置
            uq(i) = dot(B((k-1)*(Np+1)+1:k*(Np+1)),un((k-1)*(Np+1)+1:k*(Np+1))); % 计算速度场在高斯积分点处的插值值
        end
        f = 0.5 * uq.^2; % 计算速度场右端项
        u((k-1)*(Np+1)+1:k*(Np+1)) = u((k-1)*(Np+1)+1:k*(Np+1)) - dt/dx * M((k-1)*(Np+1)+1:k*(Np+1),(k-1)*(Np+1)+1:k*(Np+1)) \ (A((k-1)*(Np+1)+1:k*(Np+1),(k-1)*(Np+1)+1:k*(Np+1)) * u((k-1)*(Np+1)+1:k*(Np+1)) - B((k-1)*(Np+1)+1:k*(Np+1)) * f * dx/2); % 更新速度场
    end
    if mod(n,10) == 0 % 每隔10个时间步长输出一次结果
        plot(x,u1,'b-',x,u,'r-');
        axis([0 L -1 1]);
        xlabel('x');
        ylabel('u');
        title(['one-dimensional Burgers equation RKDG solution, t = ',num2str(n*dt)]);
        legend('Initial condition','Numerical solution');
        drawnow;
    end
end

其中，`gauss`函数用于计算高斯积分点和权重，`lagrange_basis`函数用于计算拉格朗日插值基函数，`d_basis`函数用于计算插值基函数的导数。在代码中，我们首先定义了空间和时间的区间和步长，然后定义了初始条件和有限元方法的相关参数。接着，我们循环遍历每个单元，计算基函数和测试函数的值，并计算刚度矩阵和质量矩阵。然后，我们将矩阵转换为向量，采用Runge-Kutta方法对时间进行离散化，循环遍历每个单元，计算速度场在高斯积分点处的插值值和右端项，并用矩阵求解速度场的更新方程。最后，我们输出数值解，并绘制数值解和初始条件的图像。