201711111 OKABE - 1

不連続ガラーキン有限要素法の有効性に関する基礎的検討
Basic study on the effectiveness of the discontinuous Galerkin finite
element method
岡部 つくしTsukushi Okabe（指導教員 山本亨輔）
Abstract - In this study, we compare the accuracy of numerical simulations of advection ofsquare waves under various conditions using the continuous Galerkin method and thediscontinuous Galerkin method and investigate more accurate and stable schemes.
Specifically, we compare the accuracy, stability, and computation time of the continuousGalerkin method with the Lagrange basis and the discontinuous Galerkin method with theLegendre basis, and the continuous Galerkin method with the Hermite basis and thediscontinuous Galerkin method with the Hermite basis. We will compare the accuracy,stability and computation time of the continuous Galerkin method with Hermite basis andthe discontinuous Galerkin method with Hermite basis to find out what kind of scheme ismore efficient.

１はじめにでなく、非線形システムの数値計算にも適用することができる。

数値シミュレーションは、対象とする現象の支

有限要素法の連立化プロセスは、ガラーキン法配方程式を離散的に解くことで実現される。最もが最も一般的である。これは重み付き残差法の一一般的な数値計算手法の一つは、有限要素法であ種であり、解くべき偏微分方程式の解と重み関数り、離散化と連立化によって特徴づけられる。先ず、をそれぞれ離散化する際に、同じ基底関数を用い離散化プロセスでは、支配方程式の定義領域を有る方法である。ガラーキン法は、CG（Continuous限個の微小要素に分割し、各要素内の解を既知のGalerkin ： 連 続 ガ ラ ー キ ン ） 法 お よ び DG基底関数で近似する。次に、ガラーキン法に基づき、（Discontinuous Galerkin：不連続ガラーキン）法の重み付き残差を全領域で積分することで、線形連二種類に大別でき、CG 法では、要素境界において立方程式を解く問題へと帰着させる。有限要素法

近似解の連続性が保証されるのに対し、DG 法では、は、構造力学分野だけでなく、流体や電磁場等の要素境界において近似解の不連続性が許容される。

様々な現象に適用できる。また、線形システムだけ有限要素法では、前者を用いるのが一般的だが、矩形波などの不連続面を持つ解に対しては DG 法の201711111 OKABE - 2

方が高精度である。そこで、本研究では、須藤らの義領域を[𝑥 , 𝑥 ]として積分すると、重み付き残差検討を参考に物理量が流れる様子を偏微分方程式の平均𝐼は次式のように表せる。の形で表した 1 次元線形移流方程式を対象に、CG法と DG 法の精度、安定性および計算コストを比𝐼=較する。CG 法の基底関数には 1 次 Lagrange 多項𝑤(𝑥)𝑅(𝑥, 𝑡)𝑑𝑥(2)式および Hermite 多項式を用いるのに対して、DG法 で は 、 直 交基 底 で あ る Legendre 多 項式 及 び

この時、重み関数𝑤(𝑥)は任意であり、どんな値Hermite 多項式を用いる。検討には、矩形波を採用をもとり得る。そのため、重み付き残差法は𝐼の値する。DG 法は前進差分法を採用するため、テイラが 0 となるように𝑢(𝑥, 𝑡)を定める方法であるといーガラーキン法およびルンゲクッタ法を適用し、える。
前章で述べた通り、CG 法では要素間において近その有用性も検討する。
先述の通り、DG 法は近似解の要素境界における似関数が連続であることより、領域全体の大きさ連続条件を緩和する手法であり、各要素を独立しを𝑙と仮定すると、CG 法における𝐼は次式のようにたものとして考えることができる。そのため、並列表せる。

計算が容易で、分散コンピューティングに適した手法であると言える。このことを考慮し、DG 法における数値計算では、より高精度かつ高速な計算

𝐼=

𝑤(𝑥)𝑅(𝑥, 𝑡)𝑑𝑥(3)が可能であると期待できる。
これに対して、DG 法では不連続であることが許２重み付き残差法の概要容されることより、各要素を独立したものとして重み付き残差法とは、微分方程式に数値解を代考えることができるため、DG 法における𝐼を入したときに発生する残差に、任意の重み関数を掛けることで得られる重み付き残差の定義領域全体において平均値が 0 とすることにより適切な解𝐼=

𝑤(𝑥)𝑅(𝑥, 𝑡)𝑑𝑥(4)を得ようとする手法である。具体的には以下のよと表せ、一つの要素に対する左右の離散点間が定うである。
𝑥と𝑡の値に依存するある物理量𝑢(𝑥, 𝑡)が微分方義領域となる。

程式の形で表される時、係数のベクトル𝒖と基底関数のベクトル𝑵により、𝑢(𝑥, 𝑡)は近似解として次式のように表せる。

３支配方程式の離散化
3.1 CG 法における離散化過程

𝑢(𝑥, 𝑡) ≅ 𝒖(𝑡) ∙ 𝑵(𝑥)(1)

前章で述べた重み付き残差法を利用して、今回扱う 1 次元移流線形移流方程式の離散化式を導出

これを解くべき微分方程式に代入したときに生する。未知のスカラー関数𝑢(𝑥, 𝑡)に関する 1 次元線じる残差𝑅に任意の重み関数𝑤(𝑥)を掛け合わせ、定形移流方程式は、移流速度𝑐を用いて以下のように表される。

201711111 OKABE - 3

𝜕𝑢𝜕𝑢+𝑐=0𝜕𝑥𝜕𝑡(5)

3.2 DG 法における離散化過程
前節に続き、今度は DG 法における 1 次元線形移流方程式の離散化式を導出する。式(5)に対して

この微分方程式に前章で記述した重み付き残差

重み付き残差法を適用すると、法を適用すると、
𝑤(𝑥)𝜕𝑢𝜕𝑢+𝑐𝑑𝑥 = 0𝑤(𝑥)𝜕𝑡𝜕𝑥𝜕𝑢𝜕𝑢+𝑐𝑑𝑥 = 0𝜕𝑡𝜕𝑥(11)(6)となる。ここで上式の左辺第 2 項を弱形式化し、となる。ここで、ガラーキン法に基づいて、重み関式(1), (7)を代入し、整理すると、数の近似解を𝑤(𝑥) ≅ 𝒘 ∙ 𝑵(𝑥)

(7)

𝑵𝑵 𝑑𝑥

𝜕𝒖+ 𝑐(𝑢 𝑵 − 𝑢 𝑵 )𝜕𝑡(12)−として、これを式(1)とともに式(6)に代入し、整理𝜕𝑵
𝑵 𝑑𝑥 𝒖 = 0𝑐𝜕𝑥すると、が得られる。CG 法における過程同様に、マトリク
𝜕𝒖𝑵𝑵𝑻 𝑑𝑥+ 𝑐𝜕𝑡𝜕𝑵𝑵𝑑𝑥 𝒖 = 0𝜕𝑥(8)ス部をそれぞれ𝐌、𝐂とすると、𝐌𝒖(𝑡) − 𝐂𝒖(𝑡) + 𝑐(𝑢 𝑵 − 𝑢 𝑵 ) = 0(13)が得られる。上式の左辺第 1 項のマトリクス部を𝐌、第 2 項のマトリクス部を𝐂とすると、上式は、と簡潔に書き表すことができる。また、この式を時間微分項について、前進差分法により離散化し、整𝐌𝒖(𝑡) + 𝐂𝒖(𝑡) = 0(9)理すると、と簡潔に書き表すことができる。また、この式を時間微分項について、後退差分法により離散化し、整𝒖(𝑡 + ∆𝑡) = 𝑰 + 𝑐∆𝑡𝐌 𝐂 𝒖(𝑡)∆𝑥(14)

理すると、∆𝑡+𝑐𝐌 (𝑢 𝑵 − 𝑢 𝑵 )∆𝑥𝒖(𝑡 + ∆𝑡) = [(𝐌 + ∆𝑡𝐂) 𝐌]𝒖(𝑡)(10)となる。となる。

201711111 OKABE - 4

本研究では、1～3 次それぞれの場合について数3.3 基底関数の定義
3.1 節及び 3.2 節で定義した𝐌と𝐂の中身は定義する基底関数に依存する。そこで本節では、本研究値計算を行うため、適宜、次数に合わせ、マトリクスを求める。

で用いる Lagrange 基底、Hermite 基底、Legendre 基底それぞれの場合における補間関数を紹介する。
Lagrange 基底における補間関数の各要素は以下

4 CG 法および DG 法の比較
4.1 Lagrange 基底および Legendre 基底による数値のようである。

計算例
𝑁

( )

(𝑋) =1 − 𝑋 (𝑖 = 𝑒)
𝑋 (𝑖 = 𝑒 + 1)
0 (𝑖 ≠ 𝑒, 𝑒 + 1)

(15)

まず、Lagrange 基底を用いた CG 法による計算と1 次の Legendre 基底を用いた DG 法による計算を比較する。互いに異なる基底を用いたが、どちらまた、Hermite 基底における補間関数の各要素は以下のように表される。の場合も単純な線形近似であるため、十分に比較が成り立つと考えられる。
計算条件は、計算領域を0 ≤ 𝑥 ≤ 1、移流速度𝑐 =𝑁( )1とし、要素数𝑛及び時間刻み幅∆𝑡に着眼して、同(𝑋)(𝑋 − 1) (𝑋 + 2) (𝑖 = 𝑒 = 2𝑒 − 1)⎧(𝑖 = 𝑒 = 2𝑒)(𝑋− 1) (𝑋 + 1)
⎪1
(𝑖 = 𝑒 = 2𝑒 + 1)=−(𝑋+1)(𝑋−2)
4⎨⎪ (𝑋 + 1) (𝑋 − 1) (𝑖 = 𝑒 = 2𝑒 + 2)
(𝑖 ≠ 𝑒 , 𝑒 , 𝑒 , 𝑒 )⎩0

条件の下、比較検討を行った。また、その際、𝑛(16)及び∆𝑡の値は、安定条件がより厳しい陽解法を採用した DG 法の条件に合わせた。数値計算の例を
図 4-1、図 4-2、図 4-3、図 4-4 に示す。

これらを基にマトリクスを求めるが、CG 法では離散点において連続であるため、要素ごとに得られたマトリクスの全体系への重ね合わせが必要となる。
また、Legendre 基底における補間関数は以下のようになる。

𝑁 (𝑋) = 1
𝑁 (𝑋) = 𝑋
1𝑁 (𝑋) = (3𝑋 − 1)
2
1
𝑁 (𝑋) = (5𝑋 − 3𝑋)2

図 4-1 CG 法（Lagrange 基底、後退差分）(17)

𝑛 = 100, ∆𝑡 = 0.00005

201711111 OKABE - 5

振動が小さくなるのに対し、DG 法では、矩形波が全体的に僅かに大きくなっていることが確認できる。また、図 4-3、図 4-4 においても、図 4-1、図4-2 で確認できた特徴と同じような特徴が見られる。
以上より、この条件下において CG 法では、数値振動は発生するものの、時間が経つにつれ安定していく。それに対して DG 法では、要素数に関係な
図 4-2 DG 法（Legendre 基底、前進差分)
𝑛 = 100, ∆𝑡 = 0.00005

く安定はするが、時間が経つにつれて矩形波全体が大きくなることが分かる。
次は精度及び計算にかかる時間を比較していく。
以上と同様に、𝑛及び∆𝑡の値に着目し、精度比較を行った。それらを要素数ごとにまとめたグラフを
図 4-5、図 4-6 に示す。

図 4-3 CG 法（Lagrange 基底、後退差分）
𝑛 = 500, ∆𝑡 = 0.000002

図 4-5 各要素数と格子点における誤差の二乗和

図 4-4 DG 法（Legendre 基底、前進差分)
𝑛 = 500, ∆𝑡 = 0.000002

図 4-1、図 4-2 より、CG 法に比べて、DG 法は明らかに数値振動が小さく、安定していることが分かる。また、CG 法では時間が経つにつれ、数値

図 4-6 各要素数における計算時間201711111 OKABE - 6

図 4-5 より、どちらも線形近似であるためか、ほ1.5t=0t=1.0
t=5.0

ぼ同じ精度で計算できていることが分かる。また、
1

図 4-6 より、計算速度は圧倒的に DG 法の方が遅いことが分かり、要素数が多くになるにつれそのu0.5差は大きくなっていくことが確認できる。しかし、一般的に、支配方程式は 1 次精度の前進差分法を0用いて解くことが難しいとされる。DG 法は、時間-0.5

コストは大きいものの、精度と安定性に優れ、前進的に並列計算できる点に優位性が認められると言

00.1

0.2

0.3

0.4

0.5

0.6

0.7

0.80.91x

図 4-8 DG 法（Hermite 基底、前進差分)える。また、これらの条件以外での計算を試みた結𝑛 = 50, ∆𝑡 = 0.000002

果、DG 法では時間刻み幅を小さくすればするほど、精度および安定性が高くなる。但し、後退差分 CG
1.5法と比べて、時間刻み幅を大幅に細かくしなけれ

t=0t=1.0t=5.0ば安定して計算することもできない。
1u4.2 Hermite 基底による数値計算例

0.5

次に、Hermite 基底を用いて、CG 法と DG 法の計算を行い、前節と同様に両者を比較検討してい0く。その結果を図 4-7、図 4-8、図 4-9、図 4-10、-0.5

図 4-11、図 4-12、図 4-13 に示す。

00.1

0.2

0.3

0.4

0.5

0.6

0.7

0.80.9x

図 4-9 CG 法（Hermite 基底、後退差分）
𝑛 = 250, ∆𝑡 = 0.00002

1.5
t=0t=1.0
t=5.0

u1

0.5

0-0.5

00.1

0.2

0.3

0.4

0.5

0.6

0.7

0.80.91x

図 4-7 CG 法（Hermite 基底、後退差分）
𝑛 = 50, ∆𝑡 = 0.0001
図 4-10 DG 法（Hermite 基底、前進差分)
𝑛 = 250, ∆𝑡 = 0.000008

1

201711111 OKABE - 7数に時間刻み幅を対応させて計算することとした。
0.018

CG法
DG法

0.016なるよう設定した。

0.014

誤差の二乗和また、その際、時間刻み幅は最も安定かつ高精度と図 4-7、図 4-8 より、後退差分 CG 法では時間に0.012

0.01比例してなまりが生じるのに対して、前進差分 DG

0.008

法では経過時間の影響が小さく、矩形波の形状を0.006保ったまま計算できることが確認できる。
図 4-9 で
0.004

0.002

は図 4-7 の時よりも数値振動およびなまりを抑え50100200

要素数ることができた。これは、要素数を増やしたためと

図 4-11 各要素数と格子点における誤差の二乗和考えられる。しかし、CG 法における計算では、図に示された時間刻み幅が安定に計算可能な最小刻み値であり、この値よりも細かくすると発散する

5CG法

4.54ことも分かった。また、図 4-10 では、計算の途中で明らかに発散していることが分かるが、時間刻3.5

み幅をここまで細かくしても安定した計算はでき計算時間32.5ないことが判明しており、また、この値よりも細か2くしてしまうと、PC メモリの上限の問題で計算不1.5可になることも分かっている。

10.50

図 4-11 では、異なる条件で計算を行っているが、50100

200

要素数

精度面においては、若干ではあるが、より細かい時

図 4-12 各要素数における計算時間間刻み幅で計算可能な DG 法の方がやや優れていることが分かる。図 4-12、図 4-13 からは、DG 法での計算時間の長さが前節よりも顕著に伺える。

14000

12000

5 DG 法におけるスキームの改善

計算時間10000

8000DG法

6000

5.1 テイラーガラーキン法の適用前章で課題となった時間刻みによる影響につい4000て、テイラーガラーキン法の考え方により、DG 法

2000のスキーム改善を試みる。ある時刻の解𝑢 =50

100

200

要素数

図 4-13 各要素数における計算時間𝑢(𝑖Δ𝑥, 𝑛Δ𝑡)が既知である時、時刻ステップにおける解𝑢= 𝑢(𝑖Δ𝑥, (𝑛 + 1)Δ𝑡)はテイラー展開により、前節と同様に、要素数及び時間刻み幅について、同条件下における両者の比較を試みたが、DG 法の条件を CG 法に適用すると解が発散するため、要素𝑢= 𝑢 + Δ𝑡

𝜕𝑢 Δ𝑡 𝜕 𝑢++ 𝑂(Δ𝑡 )2 𝜕𝑡
𝜕𝑡(18)

201711111 OKABE - 8と表せる。ここで、式(5)とその時間微分より、この離散化式を基に数値計算を行った結果を図5-1、図 5-2、図 5-3、図 5-4 に示す。また、ここでは前章において安定性面において優れていること

𝜕𝑢𝜕𝑢= −𝑐𝜕𝑥𝜕𝑡(19)がわかった Legendre 基底を用いることとする。

𝜕𝜕𝑢𝜕 𝜕𝑢𝜕 𝑢𝜕 𝑢=−𝑐= −𝑐=𝑐𝜕𝑡𝜕𝑥𝜕𝑡𝜕𝑥𝜕𝑥 𝜕𝑡となるため、これらを式(18)に代入し、高次の微小項を無視すると、

𝑢= 𝑢 − 𝑐Δ𝑡Δ𝑡 𝜕 𝑢
𝜕𝑢+𝑐2 𝜕𝑥𝜕𝑥(20)となる。これにガラーキン法を適用し、空間方向に

図 5-1 DG 法(1 次 Legendre 基底)ついて離散化する方法をテイラーガラーキン法と𝑛 = 500, ∆𝑡 = 0.00005いう。第 3 章と同様に、この式に DG 法を適用し、離散化するスキームを DTG（Discontinuous TaylorGalerkin：不連続テイラーガラーキン）法と呼称すuる。DTG 法を移流方程式に適用すると、式(20)より次式を得る。

𝒖(𝑡 + ∆𝑡) = 𝐌

𝐌 + 𝑐∆𝑡𝐂 − 𝑐Δ𝑡𝐊 𝒖(𝑡)2(21)+ 𝑐∆𝑡𝐌 (𝑢 𝑵 − 𝑢 𝑵 )図 5-2 DTG 法(1 次 Legendre 基底)𝑛 = 500, ∆𝑡 = 0.00005

この時、マトリクス𝐌、𝐂、𝐊は、

𝑵𝑵 𝑑𝑥

𝐂=𝜕𝑵
𝑵 𝑑𝑥
𝜕𝑥

𝐊=

𝜕𝑵 𝜕𝑵𝑑𝑥
𝜕𝑥 𝜕𝑥u𝐌=(22)

図 5-3 DG 法(3 次 Legendre 基底)である。

𝑛 = 250, ∆𝑡 = 0.00001

201711111 OKABE - 9

このように、これら 4 つの傾きを使って微分方1.5程式を近似する手法を 4 次のルンゲクッタ法といu1

い、この考え方を基に改良したスキームを RKDG(Runge-Kutta Discontinuous Galerkin：ルンゲクッタ0.5不連続ガラーキン)法と呼称する。このスキームで計算を行った結果を図 5-5、図 5-6 に示す。

0-0.500.1

0.2

0.3

0.4

0.5

0.60.70.80.91x1.5t=0t=1.0t=5.0

図 5-4 DTG 法(3 次 Legendre 基底)1u𝑛 = 250, ∆𝑡 = 0.00001

0.5

図 5-1、図 5-2、図 5-3、図 5-4 より、不連続面においては多少振動するものの、通常の DG 法で0は計算できない粗めの時間刻み幅に対して、DTG-0.5法であれば、1 次および 3 次の場合を問わず、安定

00.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.91x

図 5-6 RKDG 法(1 次 Legendre 基底)

した計算可能であることが分かり、DTG 法を用い

𝑛 = 500, ∆𝑡 = 0.0001

ることで時間刻み幅による制限を緩和することができるということを確認することができた。
1.5t=0t=1.0
t=5.0

4.2 ルンゲクッタ法の適用
1

前節に引き続き本節では 4 次のルンゲクッタ法
uの考え方を用いて、DG 法のスキーム改善を試みる。

0.5

今、ある微分方程式に対して、ある時刻𝑡 における値𝑢 及び微分係数𝑓(𝑡 , 𝑢 )が既知であるとすると、そこからΔ𝑡だけ進んだ時刻𝑡における値𝑢0は-0.5以下のように求めることができる。

00.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.91x

図 5-6 RKDG 法(3 次 Legendre 基底)
𝑘 = 𝑓(𝑡 , 𝑢 )𝑘 =𝑓 𝑡 +Δ𝑡Δ𝑡,𝑢 + 𝑘22

𝑘 =𝑓 𝑡 +Δ𝑡Δ𝑡,𝑢 + 𝑘22

𝑘 = 𝑓(𝑡 + Δ𝑡, 𝑢 + Δ𝑡𝑘 )𝑢=𝑢 +Δ𝑡(𝑘 + 2𝑘 + 2𝑘 + 𝑘 )6𝑛 = 250, ∆𝑡 = 0.0002

図 5-5、図 5-6 より、RKDG 法による計算は前節(23)で扱った DTG 法よりもさらに粗い時間刻み値での計算が可能であることが分かり、本手法は課題である時間刻みによる影響を大幅に軽減することができる優れた手法であると言える。しかし、他条件での計算を試みたところ、時間刻み値における条件は、要素数および次数の増加に比例することも201711111 OKABE - 10

分かっており、前進差分 DG 法による数値計算の参考文献

困難さを改善する余地は大いにあると言える。前[1]須藤大樹、山本亨輔：Hermite 型有限要素法

述した通り、今回は前節よりも荒い刻み値での計を用いた移流方程式と Euler 方程式の数値解算であるため、通常の DG 法による計算は発散す析、土木学会全国大会第 76 回年次学術講演ることが明らかであり、比較することができない会、2021

ため、通常の DG 法での計算結果は省略した。

[2]

松浦大志：Wavelet Taylor Galerkin 法による１次元 Euler 方程式の数値計算，筑波大学大学院システム情報工学研究科修士論文、6 まとめ
本研究では、1 次元線形移流方程式を対象に、精

2020
[3]

上野智久：天体物理学のための数値流体計度、安定性および計算時間について、CG 法と DG算アルゴリズム、東京大学大学院理学系研法の比較を行った。また、その過程で明らかとなっ究科物理学専攻修士論文、2019

た DG 法の時間刻みに関する問題に着目し、その計算方法を改良した。最後に、本研究で得られた知見を以下に示す。

1.前進差分 DG 法は、後退差分 CG 法と比べて精度および安定性には優れているが、計算に要する時間が圧倒的に長く、計算コストにおいては劣っている。

2.前進差分 DG 法では、安定した計算をする上での条件が厳しく、後退差分 CG 法に対して、大幅に時間刻み幅を細かくしなければならない。

3.通常の DG 法を改良することによって得られたDTG 法および RKDG 法のスキームを用いることによって、後退差分 DG 法による計算を困難なものとする時間刻み幅による制約を緩和させることができ、ある程度粗い分割幅においても安定した計算が可能となる。