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_jx_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的afj_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 h5 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;endend%u的坐标变换if(t(1)==j-1)for i=0:2eval(sprintf('syms u%d_j_left1',i));eval(['u' num2str(i) '= u' num2str(i) '_j_left1']);endelse if(t(1)==j+1)for i=0:2eval(sprintf('syms u%d_j_right1',i));eval(['u' num2str(i) '= u' num2str(i) '_j_right1']);endelse if(t(1)==j)for i=0:2eval(sprintf('syms u%d_j',i));eval(['u' num2str(i) '= u' num2str(i) '_j']);endendendendu=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 xx_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;endendend
else if(dif==0)if(v==q1)%因为v取q1的时候才会有±1的变化，其余情况取q0和q2都为1re=Pandingu(a,b);qv=Pandingq(re(1),re(2));elseqv=1;endend
end
end

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

matlab代码如下：

clear all;clc
syms uh_dt x j x_jleft12 x_jright12 h u0_j u1_j u2_jx_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_jx_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                2h         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;endend%u的坐标变换if(t(1)==j-1)for i=0:2eval(sprintf('syms v%d_j_left1',i));eval(['v' num2str(i) '= v' num2str(i) '_j_left1']);endelse if(t(1)==j+1)for i=0:2eval(sprintf('syms v%d_j_right1',i));eval(['v' num2str(i) '= v' num2str(i) '_j_right1']);endelse if(t(1)==j)for i=0:2eval(sprintf('syms v%d_j',i));eval(['v' num2str(i) '= v' num2str(i) '_j']);endendendendu=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_jx_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             h6 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                   h3 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                   h3 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                    h3 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=( $q^{0}$ , $q^{1}$ , $q^{2}$ )	m	l	k	m	$q^{0}$	$q^{1}$	$q^{2}$	q^3	$q^{0}_{x}$	$q^{1}_{x}$	$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值。

DG半离散格式的转化---基于matlab编写相关推荐

matlab浊音段和清音段,基于Matlab编写的语音端点检测1
wavread 基于Matlab编写的语音端点检测专业: 班级: 姓名: 指导教师: 2011年6月18日一.实验目的 1．学会MATLAB的使用,掌握MATLAB的程序设计方法: 3．掌握语音处 ...
matlab半物理仿真,一种基于Matlab的半物理仿真方法与流程
本发明属于物理建模技术领域,具体涉及一种基于Matlab的半物理仿真方法. 背景技术: 当今汽车行业,经过多年探索,业界普遍采用基于模型的控制器开发"V"模式,如图5所示.该模式可 ...
基于MATLAB编写的GNSS_SDR（GNSS软件接收机）——自学笔记（3）
今天这个跟踪函数比之前的捕获函数更有难度,我看了整整两天,才弄明白了大部分程序的含义,以下是学习记录(程序中涉及的次重点函数将在后文展示,详见目录). 目录 Tracking.m跟踪函数 ① 初始化结 ...
基于MATLAB编写的GNSS_SDR（GNSS软件接收机）——自学笔记（4）
继续来看伪距测量函数postNavigation.m函数和plotNavigation.m,其中伪距测量函数又包含了位置计算函数和伪距计算函数(leastSquarePos.m和calculatePs ...
基于MATLAB编写的GNSS_SDR（GNSS软件接收机）——自学笔记（2）
上一篇中,由于初始化中涉及到了probeData()函数,因此本文就先了解一下probeData()函数,然后是重点,介绍了捕获.捕获后的绘图函数. 目录 probeData(varargin)函数 ...
ARIMA-GRNN模型的发病率预测GUI：基于Matlab编写（ARIMA部分）
该软件首先通过经典的ARIMA模型得出初步的预测数值,生成绝对误差序列,然后输入GRNN模型得出预测的绝对误差数值,最后通过反算生成最终预测数值. 第一步.Matlab支持插件的安装 1. 双击MCR ...
ARIMA-GRNN模型的发病率预测GUI：基于Matlab编写（GRNN部分）
1.界面 1.1 ARIMA拟合值输入由ARIMA模型所取得的拟合值,如图一. 1.2实际值输入与ARIMA拟合值对应的实际值. 1.3 寻找最优光滑因子界面 1.3.1 ARIMA拟合值输出随 ...
基于Matlab的LDPC码性能研究毕业设计(含源文件)
欢迎添加微信互相交流学习哦! 项目源码:https://gitee.com/oklongmm/biye 本科毕业设计(论文) 题目 LDPC码性能研究摘要信道编码是数字通信系统的 ...
基于matlab的自动识别谱峰的程序设计,基于matlab的自动识别谱峰的程序设计毕业论文-资源下载人人文库网...
基于matlab的自动识别谱峰的程序设计毕业论文目录摘要1一绪论211几种常用寻峰方法的简单说明212小波变换413MATLAB小波分析工具箱6二小波分析基本原理721一维连续小波分析722一维离 ...
语音端点检测 matlab 论文,基于MATLAB的语音端点检测
求助,哪位高手帮忙看看以下程序全不? 基于Matlab编写的语音端点检测程序 function [x1,x2] = vad(x) %幅度归一化到[-1,1] x = double(x); x = x ...

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

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

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

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

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

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

意义

附录

DG半离散格式的转化---基于matlab编写相关推荐

最新文章

热门文章