DG半離散格式轉換微分方程組

基于下面四種類型的例子，來編寫matlab程序進行轉換：

例1、標量守恒方程的轉換

????????如下圖為標量守恒律的形式：

????????通過Galerkin法轉化為下面的變分形式：

????????其中的the global Lax-Friedrichs flux：

????????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格式)

????????

????????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格式)

????????這里假設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=(,,)	m	l	k	m				q^3				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编写的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。