DG半离散格式的转化---基于matlab编写

So it＇s me

已于 2022-10-15 22:03:09 修改

阅读量761

点赞数 2

分类专栏：数值计算文章标签： p2p linq gnu

于 2022-03-04 01:55:19 首次发布

本文链接：https://blog.csdn.net/qq_44210272/article/details/123267708

版权

数值计算专栏收录该内容

2 篇文章 2 订阅

订阅专栏

该博客详细介绍了如何使用MATLAB将不同的离散格式转换为常微分方程组，涉及标量守恒方程、通量不含对x的偏导、通量含有对x的偏导以及多个通量组合的LDG格式。通过Galerkin法和特定的基函数选择，实现了半离散格式到微分方程组的转换，为后续的时间离散和数值求解提供了便利。

摘要由CSDN通过智能技术生成

DG半离散格式转换微分方程组

基于下面四种类型的例子，来编写matlab程序进行转换：

例1、标量守恒方程的转换

如下图为标量守恒律的形式：

${{u}_{t}}+f{{(u)}_{x}}=0,\quad u(x,0)={{u}_{0}}(x)$

通过Galerkin法转化为下面的变分形式：

其中的the global Lax-Friedrichs flux：

${{\widehat{f}}_{j\pm \frac{1}{2}}}=\frac{1}{2}\left( f\left( u_{{{h}_{j\pm \frac{1}{2}}}}^{-} \right)+f\left( u_{{{h}_{j\pm \frac{1}{2}}}}^{+} \right)-{{\alpha }_{f}}\left( u_{{{h}_{j\pm \frac{1}{2}}}}^{+}-u_{{{h}_{j\pm \frac{1}{2}}}}^{-} \right) \right)$

matlab代码如下：

clear all;clc
syms uh_dt x j x_jleft12 x_jright12 h u0_j u1_j u2_j

x_jright12 = x_jleft12+h;
xj=(x_jleft12+x_jright12)/2;

%基函数的选定
% v=x.^0;%v=1
% v=2*(x-xj)/h;
v=(6*(x-xj).^2)/(h.^2)-1/2;

dv=diff(v,'x',1);

%编写基函数
q0=1;
q1=2*(x-xj)/h;
q2=(6*(x-xj).^2)/(h.^2)-1/2;

f(x)=x;
u=u0_j*q0+u1_j*q1+u2_j*q2;%u由基函数表出

re_u_dt=int(uh_dt*v.^2,x,x_jleft12,x_jright12);

a=0.5;%设置the global Lax-Friedrichs flux的a

fj_left12=0.5*(f(SpendU(-1,j-1/2,j,0))+f(SpendU(1,j-1/2,j,0))-a*(SpendU(1,j-1/2,j,0)-SpendU(-1,j-1/2,j,0)));
fj_right12=0.5*(f(SpendU(-1,j+1/2,j,0))+f(SpendU(1,j+1/2,j,0))-a*(SpendU(1,j+1/2,j,0)-SpendU(-1,j+1/2,j,0)));

re=(int(f(u)*dv,x,x_jleft12,x_jright12)+fj_left12*SolveV(1,j-1/2,v,0)-fj_right12*SolveV(-1,j+1/2,v,0))/diff(re_u_dt,'uh_dt');
re=expand(re);%u_dt 上标由v的选择决定
pretty(re)

上面求出来的结果如下：

5 u1_j   5 u0_j   5 u2_j   15 u0_j_left1   15 u1_j_left1   15 u2_j_left1   5 u0_j_right1   5 u1_j_right1
------ - ------ - ------ + ------------- + ------------- + ------------- - ------------- + -------------
   h       2 h      2 h         4 h             4 h             4 h             4 h             4 h

     5 u2_j_right1
   - -------------
          4 h

上面为d(u2)/dt的结果(程序中left表示-，right表示+)。变换v的选定，同样可以得到d(u0)/dt和d(u1)/dt，从而得到d(u)/dt=F(u)这种形式的微分方程组。

SpendU.m

function [u]=SpendU(a,b,j,d)
    %输入: u的极限属性a,a为1代表+(右极限)，a为-1代表-(左极限)；
    %      u(x)中x的下标，即x所在的单元位置；
    %      u的导数阶数d，d为1代表u有一阶导数，d为0代表u无导数。
    %输出: 由勒德让基函数表出的u，系数数值化。
    syms u0 u1 u2 h x 
    t=Pandingu(a,b);
    %基函数的计算
     if(d==1)
        q0=0;
        q1=2/h;
        q2=Pandingq(t(1),t(2))*6/h;
    else if(d==0)
             q0=1;
             q1=Pandingq(t(1),t(2));
             q2=1;
        end
     end
    %u的坐标变换
    if(t(1)==j-1)
           for i=0:2
                 eval(sprintf('syms u%d_j_left1',i));
                 eval(['u' num2str(i) '= u' num2str(i) '_j_left1']);
            end
    else if(t(1)==j+1)
            for i=0:2
                eval(sprintf('syms u%d_j_right1',i));
                eval(['u' num2str(i) '= u' num2str(i) '_j_right1']);
            end
        else if(t(1)==j)
                for i=0:2
                    eval(sprintf('syms u%d_j',i));
                    eval(['u' num2str(i) '= u' num2str(i) '_j']);
                end
            end
        end
    end
    u=u0*q0+u1*q1+u2*q2;
end

Pandingu.m

function [re] = Pandingu(a,b)
%1代表+；-1代表-。
if(a==1)
    re(1)=b+1/2;
else if(a==-1)
        re(1)=b-1/2;
    end
end
re(2)=b;
end

Pandingq.m

function [result] = Pandingq(a,b)
a_l=a-1/2;
a_r=a+1/2;
if(b==a_r)
    result=1;
else if(b==a_l)
        result=-1;
    end
end
end

SolveV.m

function[qv] = SolveV(a,b,v,dif)
%输入：基函数q的极限属性a,a为1代表+(右极限)，a为-1代表-(左极限)；
%      q(x)中x的下标，即x所在的单元位置；
%      q的导数阶数dif，dif为1代表u有一阶导数，dif为0代表u无导数。
%输出：由三个输入参数决定的q的数值。
syms h x_jleft12 x

x_jright12 = x_jleft12+h;
xj=(x_jleft12+x_jright12)/2;

q0=1;
q1=2*(x-xj)/h;
q2=(6*(x-xj).^2)/(h.^2)-1/2;

if(dif==1)
    if(v==q2)%因为v取q2的时候才会有±6/h的变化，其余情况取q0为0，q1为2/h。
        re=Pandingu(a,b);
        qv=Pandingq(re(1),re(2))*6/h;
    else if(v==q1)
            qv=2/h;
        else if(v==q0)
                qv=0;
            end
        end
    end
else if(dif==0)
        if(v==q1)%因为v取q1的时候才会有±1的变化，其余情况取q0和q2都为1
            re=Pandingu(a,b);
            qv=Pandingq(re(1),re(2));
        else
            qv=1;
        end
    end
end
end

例2、转化下面的形式为常微分方程组(通量不含对x的偏导)

matlab代码如下：

clear all;clc
syms uh_dt x j x_jleft12 x_jright12 h u0_j u1_j u2_j

x_jright12 = x_jleft12+h;
xj=(x_jleft12+x_jright12)/2;

%基函数的选定
% v=x.^0;%v=1
% v=2*(x-xj)/h;
v=(6*(x-xj).^2)/(h.^2)-1/2;
dv=diff(v,'x',1);

q0=1;
q1=2*(x-xj)/h;
q2=(6*(x-xj).^2)/(h.^2)-1/2;
f(x)=x;
u=u0_j*q0+u1_j*q1+u2_j*q2;

re_u_dt=int(uh_dt*v.^2,x,x_jleft12,x_jright12);

re=(-int(f(u)*dv,x,x_jleft12,x_jright12)+SpendU(1,j+1/2,j,0)*SolveV(-1,j+1/2,v,0)-SpendU(1,j-1/2,j,0)*SolveV(-1,j-1/2,v,0))/diff(re_u_dt,'uh_dt');

re=expand(re);%u_dt 上标由v的选择决定
pretty(re)

运行结果如下：

5 u0_j_right1   5 u1_j   5 u2_j   5 u0_j   5 u1_j_right1   5 u2_j_right1
------------- - ------ - ------ - ------ - ------------- + -------------
      h            h        h        h           h               h

上面为d(u2)/dt的结果(程序中left表示-，right表示+)。变换v的选定，同样可以得到d(u0)/dt和d(u1)/dt，从而得到d(u)/dt=F(u)这种形式的微分方程组。

例3、转化下面的形式为常微分方程组(通量含有对x的偏导ULDG格式)

${{\left( {{w}_{h}},q \right)}_{j}}-{{\left( {{u}_{h}},{{q}_{xx}} \right)}_{j}}-{{\tilde{u}}_{x}}{{q}^{-}}{{}_{j+{1}/{2}\;}}+{{\tilde{u}}_{x}}{{q}^{+}}{{}_{j-{1}/{2}\;}}+\hat{u}q_{x}^{-}{{}_{j+{1}/{2}\;}}-\hat{u}q_{x}^{+}{{}_{j-{1}/{2}\;}}=0$

$\hat{u}=u_{h}^{+},{{\tilde{u}}_{x}}=\left( {{u}_{h}} \right)_{x}^{+}$

matlab代码如下：

clear all;clc
syms uh_dt x j x_jleft12 x_jright12 h u0_j u1_j u2_j

x_jright12 = x_jleft12+h;
xj=(x_jleft12+x_jright12)/2;

% v=x.^0;%v=1
v=2*(x-xj)/h;
% v=(6*(x-xj).^2)/(h.^2)-1/2;

d2v=diff(v,'x',2);

q0=x.^0;
q1=2*(x-xj)/h;
q2=(6*(x-xj).^2)/(h.^2)-1/2;

f(x)=x;
u=u0_j*q0+u1_j*q1+u2_j*q2;

re_u_dt=int(uh_dt*v.^2,x,x_jleft12,x_jright12);

re=(int(u*d2v,x,x_jleft12,x_jright12)+SpendU(1,j+1/2,j,1)*SolveV(-1,j+1/2,v,0)-SpendU(1,j-1/2,j,1)*SolveV(1,j-1/2,v,0)-SpendU(1,j+1/2,j,0)*SolveV(-1,j+1/2,v,1)+SpendU(1,j-1/2,j,0)*SolveV(1,j-1/2,v,1))/diff(re_u_dt,'uh_dt');

re=expand(re);%u_dt 上标由v的选择决定
pretty(re)

运行结果如下：

6 u0_j   12 u2_j   6 u0_j_right1   12 u1_j_right1   24 u2_j_right1
------ - ------- - ------------- + -------------- - --------------
   2         2            2               2                2
  h         h            h               h                h

上面为d(u1)/dt的结果(程序中left表示-，right表示+)。变换v的选定，同样可以得到d(u0)/dt和d(u2)/dt，从而得到d(u)/dt=F(u)这种形式的微分方程组。

由于程序所得出的结果是化简后的最简式，若想要将半离散格式转化为矩阵格式，需要有成对的变量出现，例如出现a*u1_j，就需要有±b*u1_j±1配对，才能生成矩阵的格式，而且其系数a和b最好相同(会生成元素为1的稀疏矩阵)。此时就要对变量进行适当的加减配对。

例如：

结果只有60*u1_j，那么就对半生成(30*u1_j+30*u1_j+1)+(30*u1_j-30*u1_j+1)；

结果为20*u1_j+40*u1_j+1，那么加减配凑为(30*u1_j+30*u1_j+1)+(-10*u1_j+10*u1_j+1)。

例4、转化下面的形式为常微分方程组(含两个以上的通量组合LDG格式)

$\left( {{u}_{t}},v \right)=\widetilde{\mathcal{L}}(b(u)p,v)-{{a}_{0}}\mathcal{L}(q,v)$

$\mathcal{L}(q,v)=-\sum\limits_{j=1}^{N}{\left[ {{\left( q,{{v}_{x}} \right)}_{j}}-{{{\hat{q}}}_{j+\frac{1}{2}}}v_{j+\frac{1}{2}}^{-}+{{{\hat{q}}}_{j-\frac{1}{2}}}v_{j-\frac{1}{2}}^{+} \right]}$

$\widetilde{\mathcal{L}}(b(u)p,v)=-\sum\limits_{j=1}^{N}{\left[ {{\left( b(u)p,{{v}_{x}} \right)}_{j}}-{{(\widehat{b(u)}\hat{p})}_{j+\frac{1}{2}}}v_{j+\frac{1}{2}}^{-}+{{(\widehat{b(u)}\hat{p})}_{j-\frac{1}{2}}}v_{j-\frac{1}{2}}^{+} \right]}$

$\hat{q}={{q}^{+}},\quad \hat{p}={{p}^{+}},\quad \hat{u}={{u}^{-}},\quad \widehat{B(u)}=B\left( {{u}^{-}} \right)$

这里假设b(u)=u，在原来的基础上多出了p和q这两个变量，需要再创建这两个变量:

function [u]=SpendUAdd(a,b,j,d)
    % 创建p  
    syms v0 v1 v2 h x 
    
    t=Pandingu(a,b);
    %基函数的计算
     if(d==1)
        q0=0;
        q1=2/h;
        q2=Pandingq(t(1),t(2))*6/h;
    else if(d==0)
             q0=1;
             q1=Pandingq(t(1),t(2));
             q2=1;
        end
     end
    %u的坐标变换
    if(t(1)==j-1)
           for i=0:2
                 eval(sprintf('syms v%d_j_left1',i));
                 eval(['v' num2str(i) '= v' num2str(i) '_j_left1']);
            end
    else if(t(1)==j+1)
            for i=0:2
                eval(sprintf('syms v%d_j_right1',i));
                eval(['v' num2str(i) '= v' num2str(i) '_j_right1']);
            end
        else if(t(1)==j)
                for i=0:2
                    eval(sprintf('syms v%d_j',i));
                    eval(['v' num2str(i) '= v' num2str(i) '_j']);
                end
            end
        end
    end
    u=v0*q0+v1*q1+v2*q2;
end

代码同上(里面所有的v换成Q)，再创建q生成SpendUAdd1.m文件。主要的运行过程如下：

clear all;clc
syms uh_dt x j x_jleft12 x_jright12 h u0_j u1_j u2_j v0_j v1_j v2_j Q0_j Q1_j Q2_j

x_jright12 = x_jleft12+h;
xj=(x_jleft12+x_jright12)/2;

% v=x.^0;%v=1
v=2*(x-xj)/h;
% v=(6*(x-xj).^2)/(h.^2)-1/2;

dv=diff(v,'x',1);


q0=x.^0;
q1=2*(x-xj)/h;
q2=(6*(x-xj).^2)/(h.^2)-1/2;

a=0.5;
f(x)=x;

u=u0_j*q0+u1_j*q1+u2_j*q2;
u1=v0_j*q0+v1_j*q1+v2_j*q2;
u2=Q0_j*q0+Q1_j*q1+Q2_j*q2;


re_u_dt=int(uh_dt*v.^2,x,x_jleft12,x_jright12);

re=(-int(f(u)*u1*dv,x,x_jleft12,x_jright12)+f(SpendU(-1,j+1/2,j,0))*SpendUAdd(1,j+1/2,j,0)*SolveV(-1,j+1/2,v,0)-f(SpendU(-1,j-1/2,j,0))*SpendUAdd(1,j-1/2,j,0)*SolveV(1,j-1/2,v,0)+a*(int(u2*dv,x,x_jleft12,x_jright12)-SpendUAdd1(1,j+1/2,j,0)*SolveV(-1,j+1/2,v,0)+SpendUAdd1(1,j-1/2,j,0)*SolveV(1,j-1/2,v,0)))/diff(re_u_dt,'uh_dt');

re=expand(re);%u_dt 上标由v的选择决定
pretty(re)

运行结果如下：

3 Q0_j   3 Q1_j   3 Q2_j   3 Q0_j_right1   3 Q1_j_right1   3 Q2_j_right1   6 u0_j v0_j   2 u1_j v1_j
------ + ------ - ------ - ------------- + ------------- - ------------- - ----------- - -----------
  2 h      2 h      2 h         2 h             2 h             2 h             h             h

     6 u2_j v2_j   3 u0_j_left1 v0_j   3 u0_j_left1 v1_j   3 u0_j_left1 v2_j   3 u1_j_left1 v0_j
   - ----------- + ----------------- - ----------------- + ----------------- + -----------------
         5 h               h                   h                   h                   h

     3 u1_j_left1 v1_j   3 u1_j_left1 v2_j   3 u2_j_left1 v0_j   3 u2_j_left1 v1_j   3 u2_j_left1 v2_j
   - ----------------- + ----------------- + ----------------- - ----------------- + -----------------
             h                   h                   h                   h                   h

     3 u0_j v0_j_right1   3 u1_j v0_j_right1   3 u2_j v0_j_right1   3 u0_j v1_j_right1   3 u1_j v1_j_right1
   + ------------------ + ------------------ + ------------------ - ------------------ - ------------------
              h                    h                    h                    h                    h

     3 u2_j v1_j_right1   3 u0_j v2_j_right1   3 u1_j v2_j_right1   3 u2_j v2_j_right1
   - ------------------ + ------------------ + ------------------ + ------------------
              h                    h                    h                    h

上面为d(u1)/dt的结果(程序中left表示-，right表示+)，代码中Q代表q，v代表p。变换v的选定，同样可以得到d(u0)/dt和d(u2)/dt，从而得到d(u)/dt=F(u)这种形式的微分方程组。

意义

通过程序将DG的半离散格式快速转化为ut=F(u)的常微分方程组，解决了手动计算的失误率，以便于直接进行时间离散 Runge-kutta 的过程，大大提高了在空间离散中执行的效率。

附录

通过对基函数q的具体计算，得出了以下的数据：

计算过程的基函数数值结果
$q^{l}_{k}(x_{m})$ q=( $gif.latex?q^{0}$ , $gif.latex?q^{1}$ , $gif.latex?q^{2}$ )	m	l	k	m	$gif.latex?q^{0}$	$gif.latex?q^{1}$	$gif.latex?q^{2}$	q^3	$gif.latex?q^{3}_{x}$	$gif.latex?q^{1}_{x}$	$gif.latex?q^{2}_{x}$	q^3_x
	j-1/2	-	j-1	j-1/2	1	1	1	1	0	2/h	6/h	12/h
	j-1/2	+	j	j-1/2	1	-1	1	-1	0	2/h	-6/h	12/h
	j+1/2	-	j	j+1/2	1	1	1	1	0	2/h	6/h	12/h
	j+1/2	+	j+1	j+1/2	1	-1	1	-1	0	2/h	-6/h	12/h

此数据和上面程序的计算过程基于基函数取2个，且含非积分项含有u或q的一阶导数形式，若出现二阶及其更高阶的导数，请完善附录的数据再改变上面的代码运行。程序修改参考意见，修改基函数的数量值，请修改基函数的编写处；修改离散格式，按照具体的格式修改re；修改u或q的高阶导数，修改SpendU.m和SolveV.m文件，并分别设定相应的d和diff值。

So it＇s me

关注

2
点赞
踩
6

收藏

觉得还不错? 一键收藏
0
评论
DG半离散格式的转化---基于matlab编写

DG半离散格式转换微分方程组基于下面三个例子，来编写matlab程序进行转换：1、标量守恒方程的转换如下图为标量守恒律的形式：通过Galerkin法转化为下面的变分形式：其中的the global Lax-Friedrichs flux：matlab代码如下：clear all;clcsyms uh_dt x j x_jleft12 x_jright12 h u0_j u1_j u2_jx_...
复制链接

扫一扫

专栏目录