研0紀實：一維間斷Galerkin方法(DG)
本部分考慮簡單的一維問題：
。實際求解的算例為
並採用週期邊界條件（我導說用什麼邊界條件也是一個專門的研究方向，所以現在用什麼邊界條件不重要，能把演算法搞懂就行了）。
1。 間斷Galerkin方法（Discontinuous Galerkin Method）
假設問題的空間求解域為
![[a,b]](https://img.heatask.com/upload/website_attach/6uA56tKm1-MoT87Eo6-e76PLSvML25hgod05J+IC_FTm-2tNjzrFLNV1SSx7kNmGQ.jpeg)
，先定義網格節點：
每個單元
![[x_{i-1/2},x_{i+1/2}]](https://img.heatask.com/upload/website_attach/9CGNmCC3VrQeAZO+Ck9TGk1IFOxISdl5mpINLpAgMYmnBJZ7ouPdLr1e2FlITvvqfdNRVNIQWv6ft8sdHku6NvCDh+iWT2vXfx2O+UTj9.jpeg)
記作
，長度記為
。
DG方法與以前學的有限元方法的區別和優勢在於：
它不需要基函式是連續的，從而很大程度上提高了靈活性
。我們定義分片多項式空間
如下：
就是說
限制在每個單元上都是那個單元上的至多
次多項式。易知，它是
維的。
接下來推導方程的弱形式：設
是真解，則
![\int_0^1 {\left[ {u_t + f{{({u})}_x}} \right]{v}dx} = 0,\forall {v} \in V_h^k](https://img.heatask.com/upload/website_attach/EB1kC-APxZUbMbPXEBZ=BN1jsx=ekN=3Bo6kA5EY9g2Km6acK6x=Cvm2SJIW4f6_y+C7gGpeCN=VMqtwLciTrONzvkcRifgKzODfqU2drPzKE1HcKO1yhJEwNGHaifNf=PLfqUfdrP1-15EdlBQsBNM-60c+8Bgky+CIF-adrP1-2LtoMkinhJNUNG7IBEUolovchn6V2sEKxY4iCJn-BqWE5c=_if=+EK+e8.jpeg)
我們可以選取
的這樣一組基：它只在一個單元上非零。所以上式等價於
![\[\int_{{I_i}} {\left[ {{u_t} + f{{(u)}_x}} \right]vdx} = 0,\forall v \in V_k^h,\forall 1 \le i \le N\]](https://img.heatask.com/upload/website_attach/c17uMy25EfCz0s3TuGI9+PX_csCTO3MpfYWu990iWR6S1CikcwyAhnoTzXZD1on26U3g27Kk+mCt5k4W=-4CqaXg-Tpg639-bp6Ckool2q2k5k4a5WO4p6KUojpb6_wN6U3g27Kk2qug0xredTrvxx0hnV2xs-4w7gqveipt8QfO5kMW5Gnsp6KUojTUP17EGOqjebXPeQG_H2XU3NwuqaXy-swQt202cgqvf9RSFMOSXWbz34tQIOYD0kTUawbRb-UM3bXPeQrSu5RkcwyAhnwSc.jpeg)
利用分部積分：
![\[\int_{{I_i}} {f{{(u)}_x}v} dx = \int_{{I_i}} {vdf(u)} = - \int_{{I_i}} {f(u){v_x}dx} + \left. {f({u}){v}} \right|_{x_{i - 1/2}^ + }^{x_{i + 1/2}^ - }\]](https://img.heatask.com/upload/website_attach/GAdgf0NlFaFj1_b_Imj=Z0pvbjL-LyZw1nSgdXBTXhIx2xHTyUHTCwuDyrOEbvUeufLY7-1aIrFFuncU7T+UJzp4TzZwNO7=YBHwHQe_363RuncUa3_TC1637pekKyU3vVId-NSkVnWFHlQ17T+VI-2NUMLEJJIDufLknzoQftL9m+-ZzHkTLQ8DyiRAKYVyKUhAGTSkVnWFHlQ17T+VI-2NUMLEJJIDufLknzoQfN9xNJgy7T6iI-2N+kLjJJIDmAXwy1JC4a7dJ+vTyUyQYXG1GJtwNOnMufHj7-1aGK7jKY-ZzSkgjomDyiRwNO-quf7xfGUELP+xYgC17T+VdABDyiKEKY4qNeIc-Q5kVy+xYp-D7T7VJzp4-awhRIV2J4Y47-1aLKQxN8+DHSkZQQuDyiRwN6VqCehwy1JkVnBxYpjSa.jpeg)
這裡
在邊界上的取值就要借用隔壁的單元中的資訊了。所以這個時候，我們需要用一種叫
數值通量
的東西來代替
在邊界
處的值。
我們一般使用單調通量
，它應滿足三個條件：
1。 相容性：
。
2。 連續性：
對兩個變數都是Lipschitz連續的。
3。 單調性：
對第一個變數單調增，對第二個變數單調減。
（在理論證明的部分，我們會看到這幾個性質的具體用處）
至此，我們得到完整的弱形式：
![\[\int_{{I_i}} {{u_t}vdx - \int_{{I_i}} {f(u){v_x}dx} + {{\hat f}_{i + 1/2}}v(x_{i + 1/2}^ - ) - {{\hat f}_{i - 1/2}}v(x_{i - 1/2}^ + )} = 0,\forall v \in V_k^h\space\space(*) \]](https://img.heatask.com/upload/website_attach/8Fg8qAAIF-EyqHZ0G=kt8t_08Dn5FNoUM7s8onE3XF+Rrsm5kesSxba5VbD6M4fd1eX=mfGFI+EsMjav2cq-2T_mfxjHf4odM-jMaA9UHQjjmCSLeAgiMh65V=bvVB2+pCS4TBfsVlf=Mjal8WdYkt_5AVSvVB20pNC4TB+UHQ6yZE562cMk2T_mfxjHfZKHdCNCrUQ93aGTKCSme=+l1Zbl4-vvA6rO1eXumfG=Gk6XTLyZ2cqk8Y15AebnAeKOeab4TBfsVc6=v_55FY+lDY15V=bvVB1tx+D+7ofo9k6yZ=p5kEMt2T86AYjC+b1txAz4TB+H3a2zr=p5kEMt2T86AYjC+b1txAz4TI6o3a3xmCSAVbdzej163ZDt4WUtEKp4TIV9MF3ao9T6mWPSxblJMx3wrN_tx+T+mfptmFxa6GS5keSZLsvy9YjCWUKOddjF-ejsVcARSRi4k.jpeg)
對所有的單元
成立。
2。 數值求解
取定
的一組基為
，
，其中
是第
個單元上的一組基。對於
，數值通量可以簡單地取做迎風通量
。
在每個單元上設
（注意隨著
的變動，
都是
的函式），那麼弱形式等價於
寫成方程組即為：
其中
![\[A = \left[ {\begin{array}{*{20}{c}} {\displaystyle\int_{{I_i}} {{\phi _{i,0}}{\phi _{i,0}}dx} }& \cdots &{\displaystyle\int_{{I_i}} {{\phi _{i,k}}{\phi _{i,0}}dx} }\\ \vdots &{}& \vdots \\ {\displaystyle\int_{{I_i}} {{\phi _{i,0}}{\phi _{i,k}}dx} }& \cdots &{\displaystyle\int_{{I_i}} {{\phi _{i,k}}{\phi _{i,k}}dx} } \end{array}} \right]\]](https://img.heatask.com/upload/website_attach/v-p-FTsHQqyWvwukt_64vw1BtU_ZvB4+=5Z-D4uO7XpesbQZfMhiI2N_P6HnZTZhI0dYl47EC1G6E6LUqnLA8VNy=UbWMdDS5ul7=5+hKkoeeTlZiHZiIxCTT+HxZT0B8Sgbn-YhKkojDMvfE9H_uUd-ER1pr7KEl5MbPmIits-x16LUqnLC85dqlLHxZT0B9Jxbn-LhKk=eeGeG5jGu4p1EmQ2bsJRfAGudPmq_C1u6DMHavj64Bkjy7TBzMdCijJ9xHIdBMwpXDMvbFJLC961T7pYbSfKYuTy-F7PhRm1eeTlZfMeJvs6cDReVpCMNsT9bnHS3=hpBDMvfE9vTLuQ0OMldMdD_Bn9xGqdBLL9kpwcO58VLA0EjlLH7YBZhIUobn-YhKk=eeGeG5jGu4p1EmQ2bsJRfAGudPmq_sxJeeTpzE9vjA05zOM3cmxZhI0dwlSuElqxehr2ZfRliIxMyTruhMd4i=SpXATs5C1GfDMHaFJLC891lmULzpCIy5uWnATgys09kpUXksbyiIxCy7TBBKB3hIUobn-Y5C1u6DMvfE9H_tUr0P-ibSFC3AGlc-IdBM79Wq+LUqnLA8WJHlLZEMdDSvm97GwPhKkoeeTq0teLC96ay7T9bsZv5AGGc=wZPDxgjDMl-E9vTA05z=RbnpD2Vlud-B4fhC1Gf43US2JLC891Em1bEnc0B9J9xHIPhKk=eeTlZfMe+vU42s6HxYl3hj3dhPmq_C1u+DMvfE9H_tUr0P-ibSFC3AGlcD7dBM79WqRqwE9vjBp1EG+YbSfKhsSobn-MVlhl8sLLUreLC96ay7GJ7n6D2=J9UHIdFLq9kqeWY7.jpeg)
（太大了，全寫出來裝不下）
![\[H = - \left[ {\begin{array}{*{20}{c}} {{\phi _{i - 1,0}}\left( {{x_{i - 1/2}}} \right){\phi _{i,0}}\left( {{x_{i - 1/2}}} \right)}& \cdots &{{\phi _{i - 1,k}}\left( {{x_{i - 1/2}}} \right){\phi _{i,0}}\left( {{x_{i - 1/2}}} \right)}\\ \vdots &{}& \vdots \\ {{\phi _{i - 1,0}}\left( {{x_{i - 1/2}}} \right){\phi _{i,k}}\left( {{x_{i - 1/2}}} \right)}& \cdots &{{\phi _{i - 1,k}}\left( {{x_{i - 1/2}}} \right){\phi _{i,k}}\left( {{x_{i - 1/2}}} \right)} \end{array}} \right]\]](https://img.heatask.com/upload/website_attach/UphkTkVzJrfcBqCdaqlQCxzcDGf06mtgkRtkRtxatYB6C0XJFjII83XygCPEPDee6y3f548cStQZyO=pfAPZ3LzX0GV2-mI9ZzY7j=iqR=g6on1JF9kIVu2IpMkD7D5qPsQXAwVfvVA6on1Tei-oi1wozjh3x8Mcarnf5Vi_Yr7UUNu9hB2I8kCIpMA2tQpM-RkkV2hYYrPcHEuah+lOi1woF_L0aCKH7jxqV2PIX2-nyON6fAPX3XXQEGlqZ_y8ZzY7V2D+UxlTzo-JF9VQi1XJwEP9PD5qPsQVAX=fSiA6UAMTei-oi1woN=e2tLzcaDLQI1sqNGk6on1JF9AI8kCSzjh3xrIaUpgfIGlfvVA6UAETeihJAWNdCLIIZ_ypZzY7V2P7X2h7A6mRfLeI8kXQgC7I7i5tO83f5V3fvVA6oyqfV-QWi1WagCP9OP5qOdOCV2P7T1LJz=6JrSj8i1wKzj-4Z_C1aDQVAYi_SXlZyOzcei-oi1_JB9lfaMTHtqd_V2i+KnPc05uailPX3WhIEJA27zLeZzY7V2P7xKe6ongRfA7SV0z81gk2tL0Htqgf5V3QX2h7AKJjW8AIVmiIpMA2tQwHthqQWrsVCYQlFq2LfAP88LzZ0yP9PD5ttB3f548WSxfZAtfJFjII83KSzj-sZ_C5Zz-+TxH_YWe6ongRfA7SV0z80N22tL0Htqgf548cStQZyOzbfAPZ3LzZ0JlWZ_C5Uy3dW8nfKvjMyON6ei-Ki1wKgCPrOGzc-8OkV2hZX2-ryO=KYrlQ_V+IpMVfZ_z0UDQXBnsXNnPaHpsJV+AIVmDSzj-sZ_C5V5nf5Vi_Yr7UUNu9iykI8kCIpMA2tL0eZz-+Tw6hTkA6UAtJF9AIVmRSzjh3Tm1LkpqQV2hWX2-ryONpeihJOwEQgme2tLzcaDLQI1sqCXL6on1JF9AI83veEwQqZ_yraDQXAMsXCilzyONpWhLf_r2Iw-jDZ_C1ZzYIV2PIYrPaHqgRVMlWi1Wbzj-sZ_t0xoO_WWwfvVVTyOgKWAPZ47zZ1EPrOmvH-dgfIGGQX2-ryONpaseI8kXQgC7I7i5tPNcf5V3fvVA6on1TeihJOU0_FLA27zsHtqgQV2D+StrlyONpkUV8AtiIpMA2tL0eZz-+Tw6hTkA6oy1JFjII83C+E.jpeg)
再利用週期邊界條件，得到整體的方程組為（記
![\[C(t) = \left[ {\begin{array}{*{20}{c}} {{C_1}(t)}\\ \vdots \\ {{C_N}(t)} \end{array}} \right]\]](https://img.heatask.com/upload/website_attach/GAdgf0NlFaFj1_b_Imj=Z0pvbjL-LyZw1nSgdXBTXhIx2xHTyUHTCwvE7pJjdkVPZehwH1JC36qkHykC+Ckyj_TDyiJwN6zhmnLEe_oQfNwd0PEi7T+UI-2N7pKXJJIBNbTwy1JkVy4KuncU7T+UJzp4TJw_QozeNeI_-NoeWP+xNJxTy0TTCwODyrmAKYV0JOq8e4Kf4a7iJc-yTylTC1uDyiKxXOs=YBcwHQei363RuncUnykyj39vcaw_QBNP0ESQ7-1W36HVvc-yT37=Y3f17pWkJJZluf7ka.jpeg)
）：
![\[\left[ {\begin{array}{*{20}{l}} A&{}&{}&{}&{}&{}\\ {}&A&{}&{}&{}&{}\\ {}&{}&A&{}&{}&{}\\ {}&{}&{}& \ddots &{}&{}\\ {}&{}&{}&{}&A&{}\\ {}&{}&{}&{}&{}&A \end{array}} \right]C](https://img.heatask.com/upload/website_attach/fiU511GvNc9_Ij1gK0cs8Hoi84MUO5VLJGN5HYMjb_RwJkgQgkBRxU5SVxP-HZZ_oqZzsgxhK+Uh-_vQe5Th1T1k8Z6dGakVg_-BZ2K-O8EwzhuW1P7yeCozfHkxa3QVg_-mK-xaZ+U_-IpPwP7Xxiozfjxj7PUfgKVsKjxh6YUczIpPUQ7yeCoYVjxj6QURwbVuKgxzKcswzhuQg3+REjmx3UK_swdboqINvUKgvndYno2Sow3P1T1k3UfVswIQoqIzvUfwvnVxno9Rph7ydioz+OxdaQURUKVs6OxaZ8Ew0KmQg33Rx0JSARaxa3IVg_-BORa-ZgIwzhKQgkBRxUNl3UfZswdXoq7vvUfWvndknouOowHx1T1-3UK_thUjUQUZ05TbscUczIpPUQ7yeCoYVjxj6QURwbVuKgxzKcswzhuQg3+REjMSV2Kxa3QV17cBZ2K-ZgQw0KmQg33Rx0JSARbYswIQoqIzvUfwvnVxno9Rph7ydioz+OxdaQURUKVs6OxaZ+U_-IpPwP7Xxiozfjxj7PUfgKVsKjxh6YUcz9ktowqSrWjo3UfZGaJfdHYBZ2S-ZgQAno9RI0Tu8HJSVxTkswIRoq7t2Oxaa7sw09utowqS9G8BGOxi6QsVg_7BZxPeGLdEfIpPU0wXLmkh3UfVswdboqnyvUfWOJAwzhKQg34_1T1-3UfVt3cV17cBZ2K-ZgQw0KmQg33Rx0JSARaxa3IVg_-BORa-ZgIwzhKQ1SL11T8T3UXkt3gV17d4vUKgvndYno2Sow3P1T1k3UfVswIQoqIzvUfwvnIino2wowHT1T8T3UXkthURUKVs6OxaZ8gw0KGBow3P1T1k3UfVswIQoqIzvUfwvnIino2wowHT1T8T3UXkthURUKVs6OxaZ+U_-IpPwP7XxiSSVxP9H19_IfDBORaEvnVxdz6hJ2B=1T5y3UfZswdXoq7vvUfWvndkno9RowqS2yozfjxj7PUfgKVsKjxh6YUczIpPUQ7yeCoYVj=xakNMehG51g=-OJZznouOowHx1T1-3UXkswVWpfVsKjxh6YUczIpPUQ7yeCoYVjxj6QURwbV_ZjxhK+U_UZpWgXbREjmN4gxi65V-ebVsKHBa1PFQno2SowHT2yoPf56FH5h_oqZNKgxaazQw0K8towWT2zVSVxGxakQkb.jpeg)
再簡記為
，記
，這也就等價於
。
時間層上可以採用龍格-庫塔方法進行推進，常用的有3階格式：
%3階強保穩龍格庫塔方法
function
Qn
=
RK3
（
A，Q，t
）
k1
=
Q
+
t
*
L
（
A
，
Q
）；
k2
=
（
3
/
4
）
*
Q
+
（
1
/
4
）
*
k1
+
（
1
/
4
）
*
t
*
L
（
A
，
k1
）；
Qn
=
（
1
/
3
）
*
Q
+
（
2
/
3
）
*
k2
+
（
2
/
3
）
*
t
*
L
（
A
，
k2
）；
end
以及4階格式：
%10階4級強保穩，SSP-RK4格式
function
Qn
=
RK4
（
A，Q，t
）
Q1
=
Q
；
Q2
=
Q
；
for
i
=
1
：
5
Q1
=
Q1
+
（
t
/
6
）
*
L
（
A
，
Q1
）；
end
Q2
=
（
1
/
25
）
*
Q2
+
（
9
/
25
）
*
Q1
；
Q1
=
15
*
Q2
-
5
*
Q1
；
for
i
=
6
：
9
Q1
=
Q1
+
（
t
/
6
）
*
L
（
A
，
Q1
）；
end
Qn
=
Q2
+
（
3
/
5
）
*
Q1
+
（
t
/
10
）
*
L
（
A
，
Q1
）；
end
首先什麼是龍格庫塔方法，正常來說最自然的做法是Euler向前，即對於一個常微分方程
，我們的第一反應都是用
作為斜率，對
附近作一個線性的近似：
但是這個近似很糙，階數很低。為了獲得更高階的精度，我們把幾次上面的過程做一個加權的組合，就是龍格庫塔方法：
其中
和程式碼裡面的本質是一樣的，不過程式碼裡面使用的是所謂的“強保穩”格式，反正這樣做會大大加強穩定性，具體為什麼我暫時還不知道。
然後再解釋一下這裡的
是什麼：我們的方程是
，所以對於當前時間層的解向量
來說，
是
的一個近似，具體地位就相當於龍格庫塔方法裡面的
。
還要明確的一件事是：每個單元上的基如何選擇？一個很自然的想法是：選取
但這樣做有一個問題：我們發現
的最後一行是
當
很小時，
的最後一行會非常非常小（連MATLAB都會警告，工作精度接近奇異），所以不能選得這麼簡單。為了解決這個問題，可以對每一個基作一定的尺度變換，即除以
；另一種方法是，直接選取這個單元上的標準正交基：標準Legendre多項式，這樣一來我們的
甚至是一個對角陣，連求逆的工作都省了。
Legendre多項式的生成方式為：
所以它的導數為
標準化，再作尺度變換移動到第
個單元即為
要說明的細節就是這些了，接下來就寫正式的核心演算法吧：
首先寫一個Legendre多項式及其導數的生成程式：
%生成區間［a，b］上的標準Legendre多項式
function
［Pn，DPn］
=
Legendre
（
n，a，b
）
Q
=
@（
x
）
0。
*
x
；
DQ
=
@（
x
）
0。
*
x
；
R
=
@（
x
）
0。
*
x
+
1
；
DR
=
@（
x
）
0。
*
x
；
P
=
@（
x
）
（（
2
*
0
+
1
）
。*
x
。*
R
（
x
）
-
0。
*
Q
（
x
））
/
（
0
+
1
）；
DP
=
@（
x
）
（（
2
*
0
+
1
）
。*
（
x
。*
DR
（
x
）
+
R
（
x
））
-
0。
*
DQ
（
x
））
/
（
0
+
1
）；
if
n
==
-
1
Pn
=
Q
；
DPn
=
DQ
；
elseif
n
==
0
Pn
=
R
；
DPn
=
DR
；
elseif
n
==
1
Pn
=
P
；
DPn
=
DP
；
else
for
i
=
3
：
n
+
1
k
=
i
-
2
；
Q
=
R
；
DQ
=
DR
；
R
=
P
；
DR
=
DP
；
P
=
@（
x
）
（（
2
*
k
+
1
）
。*
x
。*
R
（
x
）
-
k
。*
Q
（
x
））
/
（
k
+
1
）；
DP
=
@（
x
）
（（
2
*
k
+
1
）
。*
（
x
。*
DR
（
x
）
+
R
（
x
））
-
k
。*
DQ
（
x
））
/
（
k
+
1
）；
end
Pn
=
P
；
DPn
=
DP
；
end
%尺度變換
c
=
（
a
+
b
）
/
2
；
h
=
（
b
-
a
）
/
2
；
Pn
=
@（
x
）
sqrt
（（
2
*
n
+
1
）
/
（
2
*
h
））
*
Pn
（（
x
-
c
）
/
h
）；
DPn
=
@（
x
）
sqrt
（（
2
*
n
+
1
）
/
（
2
*
h
^
3
））
*
DPn
（（
x
-
c
）
/
h
）；
end
初值條件（便於靈活調整，主程式裡只需統一寫成
）：
%初值條件
function
y
=
f
（
x
）
y
=
sin
（
x
）；
end
龍格庫塔的一步（一樣，便於靈活調整）：
%龍格庫塔的一步
function
DX
=
L
（
A，X
）
DX
=
-
A
*
X
；
end
接下來就是真正的主函數了！幾個值得說明的點：quadgk（f，a，b）是Gauss積分，kron是張量積。
CFL數是為了待會兒精度檢驗要用。
%間斷有限元方法求解一維雙曲守恆率：
%u_t + u_x = 0，初值條件為 u（x，0） = f（x），
%邊界條件取週期邊界，時間上使用龍格庫塔方法
function
［X，T，Q，C］
=
LegendRKDG
（
xa，xb，tb，k，kt，N，CFL
）
%輸入：空間求解域［xa，xb］，時間求解域［0，tb］，
%空間剖分段數N（時間步長自適應生成），
%分片多項式的最高次數k，基函式使用尺度變換後的標準Legendre多項式
%時間階數kt，可選用3，4階龍格庫塔方法，以及CFL數
%輸出：空間網格X，時間網格T，近似解矩陣Q，
%以及基函式的係數矩陣C
%初始化
h
=
（
xb
-
xa
）
/
N
；
X
=
xa
：
h
：
xb
；
X
=
X
‘
；
T
=
0
；
%構造剛度矩陣（由於選取了標準正交基，質量矩陣是單位陣）
G
=
zeros
（
k
+
1
，
k
+
1
）；
H
=
zeros
（
k
+
1
，
k
+
1
）；
for
i
=
1
：
k
+
1
for
j
=
1
：
k
+
1
［
fi
，
Dfi
］
=
Legendre
（
i
-
1
，
-
h
/
2
，
h
/
2
）；
fj
=
Legendre
（
j
-
1
，
-
h
/
2
，
h
/
2
）；
fjDfi
=
@（
x
）
Dfi
（
x
）
。*
fj
（
x
）；
G
（
i
，
j
）
=
-
quadgk
（
fjDfi
，
-
h
/
2
，
h
/
2
）
+
fj
（
h
/
2
）
*
fi
（
h
/
2
）；
H
（
i
，
j
）
=
fj
（
h
/
2
）
*
fi
（
-
h
/
2
）；
end
end
I
=
speye
（
N
）；
v
=
［
2
：
N
，
1
］；
W
=
I
（：，
v
）；
A
=
kron
（
I
，
G
）
-
kron
（
W
，
H
）；
%計算初值
for
p
=
1
：（
k
+
1
）
*
N
j
=
mod
（
p
-
1
，
k
+
1
）；
i
=
（
p
-
1
-
j
）
/
（
k
+
1
）
+
1
；
fp
=
Legendre
（
j
，
X
（
i
），
X
（
i
+
1
））；
g
=
@（
x
）
fp
（
x
）
。*
f
（
x
）；
F
（
p
）
=
quadgk
（
g
，
X
（
i
），
X
（
i
+
1
））；
end
%將初值裝填到C的第一列
C
=
F
’
；
C1
=
reshape
（
C
，［
k
+
1
，
N
］）；
C1
=
C1
‘
；
%開始求解
while
T
（
end
）
<
tb
t
=
CFL
*
h
^
（（
k
+
1
）
/
kt
）；
if
T
（
end
）
+
t
>
=
tb
t
=
tb
-
T
（
end
）；
end
T
=
［
T
；
T
（
end
）
+
t
］；
if
kt
==
3
C
=
［
C
，
RK3
（
A
，
C
（：，
end
），
t
）］；
elseif
kt
==
4
C
=
［
C
，
RK4
（
A
，
C
（：，
end
），
t
）］；
end
end
%利用元胞陣列儲存每一個時間層上的係數矩陣
E
=
cell
（
length
（
T
），
1
）；
for
j
=
1
：
length
（
T
）
R
=
reshape
（
C
（：，
j
），［
k
+
1
，
N
］）；
E
{
j
}
=
R
’
；
end
C
=
E
；
%這裡計算Q在每個單元左端點的值
for
j
=
1
：
length
（
T
）
S
=
C
{
j
}；
Le
=
cell
（
k
+
1
）；
for
i
=
1
：
k
+
1
Le
{
i
}
=
Legendre
（
i
-
1
，
-
h
/
2
，
h
/
2
）；
end
P
=
zeros
（
1
，
N
）；
for
i
=
1
：
N
for
q
=
1
：
k
+
1
lrd
=
Le
{
q
}；
P
（
i
）
=
P
（
i
）
+
S
（
i
，
q
）
*
lrd
（
-
h
/
2
）；
end
end
Q
（：，
j
）
=
P
；
end
Q
=
［
Q
；
Q
（
1
，：）］；
end
程式編寫好之後，輸入
［
X
，
T
，
Q
，
C
］
=
LegendRKDG
（
0
，
2
*
pi
，
4
*
pi
，
3
，
4
，
100
，
0。1
）；
來求一個時間為
![[0,4\pi]](https://img.heatask.com/upload/website_attach/EB1kC-APxZUbMbPXEBZ=BN1jsx=ekN=3Bo6kA5EY9g2Km6acK60nhJNU6JIW3Oohy+tI8.jpeg)
範圍內的解，選用三次多項式，四階龍格庫塔方法，100個單元，CFL數為0。1。求好之後，我們寫一個畫圖的程式來判斷它求的對不對。
%flash。m 顯示解隨時間變化規律
M
=
size
（
Q
，
2
）；
figure
xlabel
（
‘X’
）；
ylabel
（
‘Q’
）；
s
=
［
X
（
1
），
X
（
end
），
-
2
，
2
］；
axis
（
s
）；
［
p
，
q
］
=
size
（
Q
）；
Q2
=
zeros
（
p
，
q
）；
for
i
=
1
：
p
for
j
=
1
：
q
Q2
（
i
，
j
）
=
sin
（
X
（
i
）
-
T
（
j
））；
end
end
for
i
=
1
：
10
：
M
plot
（
X
，
Q
（：，
i
），
X
，
Q2
（：，
i
），
‘r。’
）；
axis
（
s
）；
pause
（
0。001
）；
end
在命令列輸入flash，即可看到我們求得的解：
顯示解
https：//www。zhihu。com/video/1449764962420359168
3。 理論部分
我導說這些都是套路，以後自己做的時候也得證明，所以還是整理一下。
3。1 解的存在唯一性
的弱解（弱形式的解）可能不是唯一的。但是在物理上，如果弱解滿足單元熵不等式，那麼它應該是唯一的。
單元熵不等式：
，其中
滿足
，
（它是一個原函式）。
一般來說證明單元熵不等式是比較困難的，但是可以證明：
性質1
的解
滿足下面形式的單元熵不等式：
![\[\frac{{\rm{d}}}{{{\rm{d}}t}}\int_{{I_i}} {U({u_h})} dx + {{\hat F}_{i + 1/2}} - {{\hat F}_{i - 1/2}} \le 0\]](https://img.heatask.com/upload/website_attach/Ym7zbBGvMY9NbicQ10X5HgTgG1rDF3WkbkWzZoMjyzR_cw0UXkl_ZqfQa-B0mGAQ+R3XI+fWKJVjaO0UX3HKvobSsSbB3-eF+R3XI+fWKJd-+7VwH0A_Zh+7sSbB3-eFbAqiQwxhSkUB_hlQ0Cs_ZhfQaH+7r3VDBG6UCCSEKJd-x6UbYP+NKDxG0BTu+NwzCQqiQwlqfks_9rIF5w9-vobwsSxALHA1IGGUCCT2KJd-Y6sUCXGvOnTW7_fFV1RDBG6F9f=-OOI_BzIUXkU4mm_t7hTu+A9DBGsDJW2EmIUkjIIUX3b_Zh_tsSxAMHw89AqgjWxzSw=iZ.jpeg)
其中
是平方熵，
是某個滿足相容性的數值通量。
證明
置
則因為
是
的解，對任意的
都有
。特別地，
，即
注意到
![\[\int_{{I_i}} {f({u_h}){{({u_h})}_x}dx} = \int_{{I_i}} {f({u_h})d{u_h}} \]](https://img.heatask.com/upload/website_attach/cMQ9gcc7m2XedGyJZZhxZc=vNGV6adxUfxE9yCavzJsnevboElMpCD_DMRLUhTjd6QclR=jN-NXUJ_zqdtdrJ0=4u8Cnb4ytPmILNf5DPfvnMzEoEcIpCd_DYg0nHb+HZHz5md=DGsAnNLuiZLU+Qbu5Qjzp-GutGmj6R=OEcOJbXB=fg3U+jAn+P_=e=9T1TFv5md-G6DQfJ_zLYA9wI=zFQj_+z9T1SKkZyMA4TGDeRDKoElMpCD5lP.jpeg)
及
![\[{(u_h^2)_t} = 2{u_h}{({u_h})_t}\]](https://img.heatask.com/upload/website_attach/gMEfcGanm1X5_1_nSfvng2F2y3G891BT_xGfaVYszIsEj=+82CdtHkSvH17LL_3mPHvrN-Eg6lQvP1=82JlpTwHt-g7LCHRmPHvrNHsEPeQEbi582JUHP54vH1-LL_pb_FCJRBE3S1D3=Errz.jpeg)
因此若記
，則
置
則容易驗證，這樣定義的
是滿足
的相容數值通量。現在我們可以寫
只需要證明
就好了（它的下標都是一樣的，下面省略）。由Lagrange中值定理，我們有
![\[\Theta = \tilde F({u^ + }) - \tilde F({u^ - }) - \hat f{u^ + } + \hat f{u^ - } = ({u^ + } - {u^ - })(\hat f(\xi ,\xi ) - \hat f({u^ - },{u^ + })) \ge 0\]](https://img.heatask.com/upload/website_attach/sSrUqtsERUvVll6NDuI01_mFI-OWzw3el3bUo+uUnAadDPsgJayyb9EOjmbM3S3luLVimZSYOUex9lV2tvny=j4O5q+Yd-BDBM0fGTd+MVydJskgQkL2PqWOjmbR341fvLWTmZLTNbrKlWZELXUO=+mlWfrJeMYPC6MwZDkKERzTBX7RlXrlnqWO50sylmlhC6MnZxPW6dTgufk4vwrE7AmOjmryzMYzJTqimZSYOUeoKiZA9mryb9m08LteD7HFL+VaGTdSM0edJZZ4jXUybjrO5qLYg8ozJGfz7IYOtVeoKiZE9DI0=+mRVKrSoxp1C6Mn-YPkOUex9=-kl_wvtrED8L4euMHDLLVaGTdSmReo9RZA9mryb9m08LteD7HFL+Mn-Yd+7aydJZ_0vXUctr1P8LMaK.jpeg)
其中最後一個不等號是因為
![\[\hat f(\xi ,\xi ) - \hat f({u^ - },{u^ + }) = [\hat f(\xi ,\xi ) - \hat f({u^ - },\xi )] + [\hat f({u^ - },\xi ) - \hat f({u^ - },{u^ + })]\]](https://img.heatask.com/upload/website_attach/8Fg8qAAIF-EyqHZ0G=kt8t_08Dn5FNoUM7s8onE3XF+Rrsm5kesSxba5VbDK8oeprYjF-U6GfFbsmCSAewdPdECJ4TjCWe_v2kj9ZNbyGQQAMjeG2cqkMY_L+Y254VK1e-jFZV60fH6RSRS62cMk2T_m+jjCWe_tEAzKmfpF3apSn9l93A6SEVB5VbD0Le_tEKp4TIV1LZjRsgU6EbPSxblI8Izn9ZKOk-jMZ46sVc6=v_55kEYSEeo5VbD0Le_tEaS4TIfo3a3xmCSLeWdPdDCNGj383f121eX=re6G+ZjjmCSmfSdYdT_LfEN64VKOkd22GV6GfEbaqEuq2cLg1ZWlGYjm+e_v2kjMaU6Nf-6yYJS5kedl1Zbl4TjH+UKOkdj9aU6Gf-6OZlv4k.jpeg)
根據Lagrange中值定理的要求，
在
之間。不妨設
，則根據數值通量
的單調性，顯然得到
![\[\hat f(\xi ,\xi ) \ge \hat f({u^ - },\xi ),\hat f({u^ - },\xi ) \ge \hat f({u^ - },{u^ + })\]](https://img.heatask.com/upload/website_attach/ApgbVYrikBkLzGVGOYP+VM_1RmnvrYFvNRsbxEpdLG53gv83q=r3ESauvQD7Vl5SFBjGDrRY3G0+QB_E1GR8g0CKYBjOB2KIqhDTw54J8Q9fCwyJF2REFF_GCIjtA5=Swy2nlZhJVTW3LQlMAv43qhRuKTptA5LjESzxw0R99DtL5I_3q=RJpj1uvLzt8sXIqhDfx54JVzDJj_S4FdRJdQSvT8sG+YCIxNCnlZeYi_EuQB8JgsZ5dQbvX_XAa2=X2y2noldTi_y5j_amgsU4dQSsA.jpeg)
於是
![\[\hat f(\xi ,\xi ) - \hat f({u^ + },{u^ - }) = [\hat f(\xi ,\xi ) - \hat f({u^ + },\xi )] + [\hat f({u^ + },\xi ) - \hat f({u^ + },{u^ - })]\ge 0\]](https://img.heatask.com/upload/website_attach/jT9-uR+TiIjK-k=kjKtx56DwzKaZhnO2S4d-42eKJTM8Ud-ZHyfpfPQX0B0y3i=i_qG7ffOerTqWWnHKttNeqp3NU4A7+Bpak7GUp-=3tdadVGAwTOSMsCDRlEpbizTixm=dDLkqF8DH-mHT=N8nFsDTl0A7+Bpcf-8hDLgMF849RTeS-pipiF2X0B0k5Bpcf2jbem4IurZ8I92zaN8pfPmCzz=h4lAn4IGxpNOp+L0YVGBf-tN+r4DI=4AUqougk7G7f9Oe8dZ8I3AzTOzMENBYhfM-Fnzcf8nbeSkeF84MWnHKtE8pf7vXcOcbfO0VjMXbhqzVirZ850JOyjEV4DDI10AxqYAc3HGhDLkMAI0KBJHKttN+qZDX0BAhhBpc3-8bhqzp+LL850JNzN8EENBY-FOd5.jpeg)
從而
。
□
性質2
（穩定性） 如果
具有
週期或緊支撐
的邊界，則
的解
是
穩定的，即
或
證明
由於邊界條件週期或緊支撐，將所有
上的單元熵不等式累加即得結論。□
引理
（投影定理）
稱為
到
的投影，如果
並且
。如果
是
到
的投影，那麼
![\[\left\| {Pu(x) - u(x)} \right\| \le C{h^{k + 1}}\]](https://img.heatask.com/upload/website_attach/UpPZYYpeUM6bVmA1zkE+wqt+kZIfK9jvkR5ZWEn0rcyKWuVu9uK3=ii7FuveTSTWZzQ7euXugMibP57M_4WMWttmA3mfZLgI7NAdeubHzuV7bA4M4aX+U90ae4_JYgxJaDhahbR5gZnK=9yJ_4_uagy8W9mI6QXS7yhcirjm1riZPA4M5nH2F.jpeg)
這裡
依賴於
，但和
無關。
性質3
（收斂階） 對於線性問題
，如果其解
是光滑的，那麼其收斂階估計為
其中
是弱形式
的解，
依賴於
，但和
無關。
證明
沿用之前的記號，因為
是
的解，故
由於
是真解，故顯然也有
（這是因為：弱形式對於真解一定是對的）
由
的線性性，有
記
，則將
取做
，得：
對左邊，此前已經證明
![\[{B_i}({e_h},{e_h}) = \frac{1}{2}\frac{{\rm{d}}}{{{\rm{d}}t}}\int_{{I_i}} {{{({e_h})}^2}} dx + {{\hat F}_{i + 1/2}} - {{\hat F}_{i - 1/2}} + {\Theta _{i - 1/2}}\]](https://img.heatask.com/upload/website_attach/_JS=ItuMMvGofv6vIv=GcsVL+7IBJ5hWRgA=f+4ByoUm+BsDyb+OPa70_KXhQ5gF2KRy1VUyO6nmZuTDy=FOGd00_KXNQ5=F2KRy1Ve2KCzH7iZ2TDjRaf00_KWaXw-BwBej1xbES3V7SuIzb++OPt70_KWNvk5hjAwEUvWyO6-mP4kDy=FOPt70_KWNv3kF2NOT+YbERI-mP4kDy=As6YI2rBwPXwhA_gu=cTbERlVoSYO+dYOlT_VrV2mNv3kF2KJyg=8ymEjmP4rDaDtOPt-04DiNv3AF2NwToZogKCIH7uV5nTORSwD0_KWNv3kF2NOHdIW2SlVoTN1Dy=kG7SV=k=maXwV1e=wEwjbES3DL7iZ4TxOlSwVpk7ETwPM1wBegcTbERIYcmwBk7ZkuGwVrV2z=NPLFeyJ2oZojKCZGX2-D+sH0ZSVrk+iHJULWwBKs1xbES3VoTyZ2TTOtThFHj.jpeg)
其中
待會兒累加的時候可以消去，不用擔心，而
。所以這邊沒什麼大問題。而對於右邊，把它具體寫出來就是
由投影的第一條性質知
，由第二條性質知
，從而上式簡化為
![\[ - {B_i}({\varepsilon _h},{e_h}) = - \int_{{I_i}} {{{({\varepsilon _h})}_t}{e_h}} dx \le \frac{1}{2}\left( {\int_{{I_i}} {{{({{({\varepsilon _h})}_t})}^2}} dx + \int_{{I_i}} {{{({e_h})}^2}} dx} \right)\]](https://img.heatask.com/upload/website_attach/3YrJjiXlV2YJ4evIyguwN_soyRPpcRELuEbJiGVTMJtr5oRn88co99hZaqwr1j_b0bHL1_5kWs3rbbDn880eB0i9i_cV_ik-DagXUqzW_DS8tss9BmaU_rPMTrt6w+t+prVacpKk9FGV2TrhiqjL_rP0lGPVUryXlusmwf=Q4uEJ72Me3KAHLcs6t68+lSCJ5W0NjQnUOIJXGd6n8_=oWnl9HjoQQctRprM7lgFuNGEJn9Me4bqFzcI9H9cohqJVdsffjgq6_DiGyvMe4bAf4nh9Hjnr199w2E4JUq_Rj2EJ72Mc3Q8JyZw9Hjtr1j_=0bHL1_5k9fsXtss9iqjL_rP0TrtKUryblus1Uq_R_DiGts0oPgMHC0q5gHLn_6Kv3rM7Db=eLlEJnqgbiqjq_rhnTr59UrEYg0M7Db=Q5GoC59VnVMTk_rDLgHTfz+tRo0M7lMnuNlEJn9Me4bpo9jh9Hjtr1jiVg26L1_jkAI3rbbqnV7No9j39H982UryXlusmwQuR_Di9=GMc3R-wCRHITrt6UrEWlu0m1.jpeg)
將所有單元上的不等式累加起來，得到：
![\[\frac{{\rm{d}}}{{{\rm{d}}t}}\int_0^1 {{{({e_h})}^2}} dx \le \int_{{I_i}} {{{({{({\varepsilon _h})}_t})}^2}} dx + \int_{{I_i}} {{{({e_h})}^2}} dx \le C{h^{2k + 2}} + \int_0^1 {{{({e_h})}^2}} dx\]](https://img.heatask.com/upload/website_attach/bJ=4Q5dVMv1ocoanILQp+5Tnf0GtbIykegL4_cbqyokmdGG8yT0lhXf8iTzm_cGQKB_ywVQ2KCcGbIG8yGnkpVbLqRIE-DRFKB_ywVQ2KC63ICx9fL4lhu+NqRIE-DRFe2FvEgBnS3Y7jntn+DmwpVxM3QDl0fUDOKSA2u-MmESmC6+8a0jlhu_83iYl0fDDONFiwVQGKC6HJnj6nOV1wsU8rOV1=ILveHYA2u-MO6qkWnw8yG-lhu_YqRIC-DRNKB_ywV-oKC63ICbH7R74pVbHqRZDjcGzahm__cNrclBUZ3TK-AVyikTKkDG0-DRFKBStwVQGKCc0LlTK-AV3xgvN7ADlcCU8KBNB_cnKYvYojlTKTwAf+kTKkAV3lGrDOKSA2u-MO6pmLEy8yGnldsy8iG-lcsPDOKtA2XBuKC6HICbLn0U4qjTIx0XlUfeDOKT8wVPHKC63LIuY7RnCqQf8iG-l0fD8KBTyxjB6RnNrc0qB7RZF2ku8iG7l0fUDOKSAeUuMO6qmWny8yG-l2Uw8iG-l0FCzKB_CwVQG6nkpICx97RZEh.jpeg)
設
![\[\varphi (t) = \int_0^1 {{{({e_h}(t))}^2}} {\rm{d}}x\]](https://img.heatask.com/upload/website_attach/J3ovuu4u3QFtujIws2SD0v8SfQ1AOtHM7aYv4d2_Q=ICvk5XeW2UaZtxi8isGeDLJ2MHDWqh7h7jenEXh5wFsXmkIGDumRbW_V8aEQZgk-7toIWTpSXbbi8jxgc70hLSVy8QfukT4wpCgKZXe5wUaZ8dqq6OsWvY8i8=p-Zyk=4AloJd3zXi7i8jkiSCaYhW_VnCe.jpeg)
，則上式成為
由Gronwall不等式，可得
![\[\varphi (t) \le C{h^{2k + 2}}{e^t}\]](https://img.heatask.com/upload/website_attach/osgE6tumoT7FXAo18pC5IsZAI97y29jfrfsEJ+4CFEzzYhluT_o_iaqPjilFfTYeovS5rZ8AXj_-EaOuT_lg+XCy8gWRItgEaXjBkV93dT_-zVOBmhWQqY7z8gWR+zgEa+zh_=WzDbozDdX5U.jpeg)
，從而
![\[\varphi (0) \le C{h^{2k + 2}}\]](https://img.heatask.com/upload/website_attach/+lPuU5=b2tHDkMVG7=kLw4eNTEvm6FSiKJ5uTcsmjqwPV583iesB=WMgtcTpTKbh+=jGXVmL35WS=g13ieSEUxzhZar5UR_qxahtCUGIrtWSON1EqcqNakriZaJ0ZbS3T.jpeg)
，即
因此即得
![\[\left\| {u( \cdot ,0) - {u_h}( \cdot ,0)} \right\| \le \left\| {{e_h}(0)} \right\| + \left\| {{\varepsilon _h}(0)} \right\| \le C{h^{k + 1}}\]](https://img.heatask.com/upload/website_attach/tx+kO-l3wqlcmaEaHRGgEk+ztCtelz62DsUkC5tg8Xe6GKxfVHkco4Gj7+UdqrpUC1p+iG9RhqJcej2fA8miy26QtCccLVgcRuA4iGGzhqCU24F7G3wzy29ROeGZ1Z+wLWVuNkbtg117Fp2qWMvk0ngj7+Urrr9zD76V3n+P1qi6ku7c9M+co4Uh=7pV0K6dC1G+_n+P1LJcekuCLIvk0n+VnjS+DfKcKkfQiG6RCUQR9WmfVHkcoVlsOeHn1Z+wLWnfBybQq3E6kXTQ3mOJy29UOe6ILCRnstxrC5JhBL_z9p2X+IvnpISjTaK+wPbiC1p+OS_5Axd6kuTfV-kiy26QDE+RcZ+kLUaf80+Qqk1U392qfq+2y29ROe9J0K6dC1pIr.jpeg)
（這裡不是同一個
，但總之是個常數）
再結合穩定性理論，即有
![\[\left\| {u( \cdot ,t) - {u_h}( \cdot ,t)} \right\| \le \left\| {{e_h}(t)} \right\| + \left\| {{\varepsilon _h}(t)} \right\| \le C{h^{k + 1}},\forall t\ge0\]](https://img.heatask.com/upload/website_attach/fiU511GvNc9_Ij1gK0cs8Hoi84MUO5VLJGN5HYMjb_RwJkgQgkBRxU5SVxP-HZZ_oqZCvUfUscU_-DpQ1Sb=1T8T84TT3PsV1Mz5vUK_sc0Ano2wJYFv1T1-3UKQthUjUQzZ05SEvnJxIZpWhU7yeCSSVxPd0ZdAJbVuKgxhKcswzqZieUQRxUYj9YbhswVWoqICsgxhK+U_-y9GfP7yeCoYW1SxPSkVg_-mvUXU1PllegKQgkBRx0Nl3U6ZthUjUQXBGpS-ZpMwzhOtowHx1T1k3UXk30FfeHPjHYVQ0+t4eZpPwP7XkHJSARWxa3QNoqZC1HEtHzQwzqOQg3B=1T8TM0xB6hURUgjBZxx-ZgJAodpWUQQ71T1-3UfVswJWoqZCGHMamLxTojKQgkWurTVSVxGxakQkb.jpeg)
□
4。 Numerical Result（數值結果，驗證收斂階）
我們在主演算法的程式裡已經寫了
這樣當空間步長變化時，時間步長自動減小，並且不會影響整體的收斂階。所以只要調空間的單元數就行了。
下面編寫測試指令碼，先寫一個計算各種誤差和Order的函式：
%誤差分析，對解和真解進行比較
function
［L1，L2，C，kL1，kL2，kC，meshsize］
=
analysisDG
（
xa，xb，tb，K，tk，N0，CFL，n，solve
）
%說明：
%輸入：X，T是空間和時間的網格剖分，Q為數值解矩陣，初始剖分段數N0，K0，
%加細次數n，用於求解的函式solve。
%輸出：L1模，L2模，L∞模隨時間變化的向量L1，L2，C，
%以及差比向量kL1，kL2，kC，用於驗證收斂階。
for
k
=
1
：
n
［
X
，
T
，
Q1
］
=
solve
（
xa
，
xb
，
tb
，
K
，
tk
，
N0
*
2
^
（
k
-
1
），
CFL
）；
［
p
，
q
］
=
size
（
Q1
）；
meshsize
（
k
，：）
=
［
p
-
1
，
q
-
1
］；
Q2
=
zeros
（
p
，
q
）；
for
i
=
1
：
p
for
j
=
1
：
q
Q2
（
i
，
j
）
=
sin
（
X
（
i
）
-
T
（
j
））；
end
end
K0
=
length
（
T
）；
%計算L1模
L1
（
k
）
=
sum
（
sum
（
abs
（
Q1
-
Q2
）））
/
（
N0
*
2
^
（
k
-
1
）
*
K0
）；
%計算L2模
L2
（
k
）
=
sqrt
（
sum
（
sum
（
abs
（（
Q1
-
Q2
）
。^
2
））））
/
sqrt
（
N0
*
2
^
（
k
-
1
）
*
K0
）；
%計算L∞模
C
（
k
）
=
max
（
max
（
abs
（
Q1
-
Q2
）））；
L1
=
L1
‘
；
L2
=
L2
’
；
C
=
C
‘
；
end
%計算階
kL1
=
L1
（
1
：
end
-
1
）
。/
L1
（
2
：
end
）；
kL2
=
L2
（
1
：
end
-
1
）
。/
L2
（
2
：
end
）；
kC
=
C
（
1
：
end
-
1
）
。/
C
（
2
：
end
）；
kL1
=
log
（
kL1
）
/
log
（
2
）；
kL2
=
log
（
kL2
）
/
log
（
2
）；
kC
=
log
（
kC
）
/
log
（
2
）；
end
然後就是測試指令碼：
%testDG。m，測試DG演算法的精度
N
=
10
；
N0
=
1
；
CFL
=
0。2
；
［
L1
，
L2
，
C
，
kL1
，
kL2
，
kC
，
meshsize
］
=
analysisDG
（
0
，
2
*
pi
，
2
*
pi
，
3
，
4
，
N0
，
CFL
，
N
，@
LegendRKDG
）；
n
=
1
：
1
：
N
；
plot
（
n
，
L1
，
’-b‘
，
n
，
L2
，
’-r‘
，
n
，
C
，
’-k‘
）
legend
（
’L1 Norm‘
，
’L2 Norm‘
，
’L∞ Norm‘
）
xlabel
（
’n‘
）
ylabel
（
’Error‘
）
title
（
’Error - n‘
）
m
=
1
：
1
：
N
-
1
；
plot
（
m
，
kL1
，
’-b‘
，
m
，
kL2
，
’-r‘
，
m
，
kC
，
’-k‘
）
legend
（
’L1 Norm‘
，
’L2 Norm‘
，
’L∞ Norm‘
）
xlabel
（
’n‘
）
title
（
’EOC‘
）
（注意這裡CFL數要自己調一下，不然可能結果很差）
最後畫出來一些類似於這樣的圖：
這是使用3次多項式，4階龍格庫塔方法得到的結果，數值結果說明收斂階是4階，很好地驗證了上面理論部分的正確性。
最後，看論文裡面好像最後都是做成這種表格的形式，我也試著弄了一下：
3次多項式，4階龍格庫塔方法的詳細結果
（後來發現網格那裡只要寫空間的size就行了，但也懶得改了，就這樣吧）
Reference： Chi-Wang Shu， Discontinuous Galerkin Methods： General Approach and Stability
DGpaper。pdf
710。9K
·
百度網盤
DG文獻筆記。pdf
99。6K
·
百度網盤
