Freefem++ による楕円型偏微分方程式の数値計算

Freefem++ による楕円型偏微分方程式の
数値計算
齊藤宣一 ∗
東京大学大学院数理科学研究科
2015 年 5 月 29 日
概要
有限要素法は，その抜群の汎用性の高さから，数ある偏微分方程式の数値解法の
中でも，最も強力なものの一つである．しかし，そのアルゴリズムを正確にプログ
ラミングすることには，（単調な作業である一方で）それなりの経験が必要になる．
例えば，O. Pironneau 教授は，彼の著書 “Finite Element Methods for Fluids
(Wiley, 1989)” のなかで次のように述べている：
Numerical analysis is somewhat dry if it is taught without testing the
methods. Unfortunately experience shows that a simple finite element
solution of a Laplace equation with the P 1 conforming element requires
at least 20 hours of programming time; so it is diﬃcult to reach the more
interesting applications discussed in this book in the time allotted to a
Master course. [page 197]
このような “非数学的” な障壁を超えるために，フリーソフトウエア Freefem++
http://www.freefem.org/ff++/
あるいは，その GUI 版である Freefem++-cs
http://www.ann.jussieu.fr/~lehyaric/ffcs/index.htm
の利用は有用である．Freefem++ を用いれば，プログラミングの知識はあまりな
くても，偏微分方程式の弱形式を理解していれば，有限要素法による数値計算が体
験できるのである．
この文書では，Freefem++ を用いた，楕円型偏微分方程式の数値計算の初歩を
解説し，応用として，半線形楕円型偏微分方程式の数値計算を述べる．
この文書で使用したプログラムは，
http://www.infsup.jp/saito/materials/examples.zip
により利用できる．
∗
norikazu[AT]ms.u-tokyo.ac.jp
1
目次
1
Freefem++ と gnuplot の準備
1
2
はじめの例題
4
3
gnuplot による可視化
10
4
領域と三角形分割
15
5
例題集
18
6
有限要素解の収束
22
7
非線形問題への応用：Nehari の変分原理
25
8
非線形問題への応用：スケーリング反復法
30
参考文献
35
1 Freefem++ と gnuplot の準備
1. Freefem++ のウェッブページ
http://www.freefem.org/ff++/index.htm
から自分の使用しているシステムに合わせて Freefem++ の最新バージョンをダ
ウンロードしてインストール（ファイルを展開するのみ）する．
2. なお，現在では，GUI 版の Freefem++-cs
http://www.ann.jussieu.fr/~lehyaric/ffcs/index.htm
が公開されている．この文章では，この Freefem++-cs を利用する．
3. 起動すると図 1.1 のようになる．主に利用するのは，図 1.2 に示した，「保存」，
「実行 (run)」，「停止」の 3 つのボタンである．
4. 次節で詳しく説明するが，左上のウインドウにプログラムを書き (図 1.3)，その上
にある「実行」ボタンを押すとそのプログラムが実行され，右のウインドウに結
果が表示される (図 1.4)．エラーなどがある場合には，下のウインドウに行番号
などが表示される．
5. あとで gnuplot も使うので，これも自分の使用しているシステムに合わせて用意
しておく．google などで調べること．いくらでも出てきます．
6. 最近，Freefem++ の参考書が出版されました．
[5] 大塚厚二，高石武史：有限要素法で学ぶ現象と数理—Freefem++ 数理思
考プログラミング—，共立出版，2014 年
1
図 1.1 Freefem++
保存
実行
停止
図 1.2
Freefem++ で主に利用するボタン
2
図 1.3
左上のウインドウにプログラムを書く
図 1.4 「実行」ボタンを押すとそのプログラムが実行され，右のウインドウに結果が表示される
3
2 はじめの例題
はじめに，有限要素法の用語を確認しておく．多角形領域 Ω ⊂ R2 に対して，次の 3 つ
の条件が成り立つように，Ω を小さな (閉) 三角形 T の集合 T に分割する（図 2.1–2.3）:
1. すべての三角形をあわせると Ω に一致する
2. 任意の 2 つの三角形はその内部を共有しない
3. 任意の三角形の任意の頂点は，他の三角形の頂点と一致しているか，あるいは単
独で Ω の角をなすかのどちらか（すなわち，ある三角形の辺上に他の三角形の頂
点があることは無い）．
三角形の小ささ（分割の細かさ）を表すパラメータ
(hT = T の外接円の直径) .
h = max hT ,
T ∈T
を導入して，これを T = Th と明示する．Th を Ω の三角形分割 (triangulation) と
言う．また，Th を構成する三角形を要素 (element)，各要素の頂点を節点 (node) と
呼ぶ．
図 2.1 Ω = (0, 1) × (0, 1) の三角形分割の例．規則的な格子を用いている．節点数
は左から，81，289，1089 であり，要素数は 128，512，2048 となっている．
Ω ⊂ R2 が十分に滑らかな境界を持つような有界領域の場合は，次のように三角形分
割 Th を導入する．まず，境界 ∂Ω 上に点 P1 , . . . , PM を配置して，隣り合う 2 つの点の
距離の最大値を h とする．そして，これらを結んでできる多角形 Ωh を考える．あとは，
上と同様に Ωh の三角形分割 Th を作れば良い（図 2.4，2.5）．
まとめると，
∪
Ωh =
T
T ∈Th
と書ける．そして，
{
=Ω
Ωh
̸= Ω
(Ω: 多角形領域)
(Ω: 滑らかな領域)
4
図 2.2 Ω = (0, 1) × (0, 1) の三角形分割の例．不規則的な格子を用いている．節点
数は左から，95，333，1267 であり，要素数は 156，600，2404 となっている．
図 2.3 多角形領域の三角形分割の例
図 2.4
楕円 (の内部) の三角形分割の例．領域が多角形でなくても，はじめに領域を
多角形で近似してから，三角形分割を作れば良い．
図 2.5 区分的に滑らかな領域の三角形分割の例
5
である．
Th 上の連続区分一次 (P1) 要素とは，
Xh = {vh ∈ C(Ωh ) | vh は各 T ∈ Th 上で一次多項式 },
で定義される関数空間（あるいはその元のこと）である．境界条件も考慮した
Vh = Xh ∩ H01 (Ωh ) = {vh ∈ Xh | vh |Γ = 0}
を連続区分一次要素と言うこともある．このとき，Xh は H 1 (Ω) の，Vh は H01 (Ω) の有
限次元部分空間となり，したがって，閉部分空間である．
Freefem++ の使い方を説明するために，次のような簡単な例題を考える．まず，
Ω = {(x, y) | x2 /4 + y 2 < 1},
Γ = ∂Ω = {(x, y) | x2 /4 + y 2 = 1}
(2.1)
とする．そして，Poisson 方程式の境界値問題
− ∆u = 1 − x2 y in Ω,
u = 0 on Γ
(2.2)
を考える．f (x, y) = 1 − x2 y とおくと，弱形式は，
∫∫ (
u∈V =
H01 (Ω),
Ω
∂u ∂v
∂u ∂v
+
∂x ∂x ∂y ∂y
)
∫∫
dxdy =
f v dxdy
(∀v ∈ V )
Ω
となる．
一方で，有限要素近似は，
∫∫
(
uh ∈ V h ,
Ωh
∂uh ∂vh
∂uh ∂vh
+
∂x ∂x
∂y ∂y
)
∫∫
dxdy =
f vh dxdy
(∀vh ∈ Vh )
Ωh
となる．
この問題を Freefem++ で解くプログラムが example1.edp (リスト 2.1) である．*1
リスト 2.1
example1.edp
//example1.edp
func g = 0;
func f = 1 - (x^2)*y;
int n = 40;
border G1(t = 0, 2*pi) {x = 2*cos(t); y = sin(t);}
mesh Th = buildmesh(G1(n));
fespace Vh(Th, P1);
Vh u,v;
solve poisson(u, v) =
int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) - int2d(Th)(f*v) + on(G1, u = g);
plot(u,ps="example1.eps");
6
図 2.6
example1.eps
example1.edp 実行すると，図 2.6 が出てくる．example1.edp を細かく見ていこう．
// example1.edp
func g = 0;
func f = 1 - (x^2)*y;
int n = 40;
//で始まる行はコメント行であり，計算の実行に関係ない．func で定義したものは x, y
の関数となる．g は境界条件を記述する関数，f は f (x, y) を意味する．int で定義した
変数は，整数として認識される．n はメッシュの細かさを制御するパラメータで，この値
を大きくするほどメッシュは細かくなる．
border G1(t = 0, 2*pi) {x = 2*cos(t); y = sin(t);}
mesh Th = buildmesh(G1(n));
Γ = ∂Ω を，
x = 2 cos t,
*1
y = sin t
(0 ≤ t ≤ 2π)
edp = équations aux dérivées partielles(偏微分方程式)
7
とパラメータ表示して，これを G1 として定義している．この場合，曲線の向きを左向
きにしなければならない．そうすることで，その内部が自動的に Ω として認識される．
そして，Ω の三角形分割を Th としている．∂Ω 上に n 個の節点を配置して，さらに内部
はそれに合わせて三角形要素が生成される．
fespace Vh(Th, P1);
Vh u,v;
Vh は Th 上で定義された区分一次 (P1) 要素である．関数 u, v を P1 要素，すなわち，
u, v∈ Vh としている．fespace = finite element space の意味である．
solve poisson(u, v) =
int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) - int2d(Th)(f*v) + on(G1, u = g);
ここでは，
∫∫
poisson (uh , vh ) =
Ωh
(
∂uh ∂vh
∂uh ∂vh
+
∂x ∂x
∂y ∂y
)
∫∫
dxdy −
f vh dxdy
Ωh
と定義して，
∀vh ∈ Vh
poisson (uh , vh ) = 0,
を解いて，uh ∈ Vh を求める．poisson(u, v) の，はじめの u が未知関数で，後ろの v
が試験関数である．
∂uh
dx(u) ←
,
∂x
∫∫
int2d(Th)(...) ←
... dxdy
Ωh
の意味である．最後の + on(G1, u = g) で境界条件 u|∂Ω = 0 を指定している．
plot(u, ps="example1.eps");
u を等高線表示して，その結果をファイル example1.eps に出力する．
問題 2.1. example1.edp を作成し，実行せよ．さらに，それを修正することで，次を
行え：
1. n の値を変えて実行するとどうなるか．
8
2. Ω = {(x, y) | x2 + y 2 = 4} に変更して，計算せよ．
3. f (x, y) = x − y 2 に変更して，計算せよ．
4. 最後の行を次のように変更するとどうなるか？
plot(u, Th, ps="example1.eps");
9
3 gnuplot による可視化
前節に引き続き，Poisson 方程式の境界値問題 (2.1)，(2.2) を考える．本節では，計算
結果の数値をファイルに書き出して，それを gnuplot で 3D 表示する．また，その結果
を eps 形式ファイルに保存する．そのために，次の example2.edp (リスト 3.1) を使う．
リスト 3.1
example2.edp
// example2.edp
func g = 0;
func f = 1-(x^2)*y;
int n = 80;
border G1(t = 0, 2*pi) {x = 2*cos(t); y = sin(t);}
mesh Th = buildmesh(G1(n));
fespace Vh(Th, P1);
Vh u,v;
solve poisson(u, v) =
int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) - int2d(Th)(f*v) + on(G1, u = g);
// data file for gnuplot
ofstream ff("example2.plt");
for(int i = 0; i < Th.nt ; i++){
for(int j = 0; j < 3; j++){
ff << Th[i][j].x << " "<< Th[i][j].y << " " << u[][Vh(i,j)] << endl;
}
ff << Th[i][0].x << " " << Th[i][0].y << " " << u[][Vh(i,0)]
<< endl << endl << endl;
}
// mesh
plot(Th, wait=true, ps="example2mh.eps");
変わっているのは次の部分のみである．
// data file for gnuplot
ofstream ff("example2.plt");
for(int i = 0; i < Th.nt ; i++){
for(int j = 0; j < 3; j++){
ff << Th[i][j].x << " "<< Th[i][j].y << " " << u[][Vh(i,j)] << endl;
}
ff << Th[i][0].x << " " << Th[i][0].y << " " << u[][Vh(i,0)]
<< endl << endl << endl;
10
}
まず，ofstream ff("FILENAME") でデータを出力するファイル FILENAME を指定する．
名前は任意で良いが，example2.plt とした (plt = plot の意味を込めている)．プログ
ラム内では，以後，ff が example2.plt を意味する (ff も，f1 でも ff0 でも任意で良
い)．次に，
ff << VAR;
変数 VAR の値を ff に出力
ff << "text";
文字列 text を ff に出力
ff << VAR << " " ;
変数 VAR の値を ff に出力した後，空白を入れる
ff << endl;
改行を入れる
を意味している．使用されている変数（VAR）は，次の通り：
Th.nt
Th の全要素数 NE
Th[i][j].x
i 番目の要素の j 番目の節点の x 座標
Th[i][j].y
i 番目の要素の j 番目の節点の y 座標
u[][Vh[i][j]]
i 番目の要素の j 番目の節点における uh の値
（ただし，i と j は次の範囲を動く：0 ≤ i ≤ NE − 1, j = 0, 1, 2）
example2.plt は，図 3.1（左）のように，3 つの数字が並ぶ行が 4 行続き，2 行分の
空白の後に，再び 3 つの数字が並ぶ行が 4 行ある，というファイルになっている．各ブ
ロックについて，1 行目と 4 行目はまったく同じ数字になっていることに注意せよ．ま
た，例えば，2.7483e-032= 2.7483 × 10−32 の意味である．これらの数字の意味は，図
(i)
(i)
3.1（右）の通りである．ただし，第 i 番目の要素の第 j 節点の座標を (xj , yj )，そこ
(i)
での uh の値を uj とする．
gnuplot のコマンド splot は，5 個の点 {Pi = (xi , yi , zi )}4i=0 が，図 3.2（左）のよ
うに与えられた際に，これらを順に線分で結んでいく．空白が 2 行ある場合には，その
前後の点を結ばない．したがって，この例の場合，(x2,y2,y2) と (x3,y3,y3) との間
には線は描かれない．結果は，図 3.2（右）のようになる．
実際に，gnuplot を使用してみよう．gnuplot を起動し，作業中のフォルダ (exam-
ple2.edp を保存したフォルダ) に移動し，gnuplot 1 のように入力すると，図 3.3（左）が
出てくる．
gnuplot 1:
gnuplot> splot "example2.plt" with lines
11
y0
(0)
y1
(0)
y2
(0)
y0
(0)
u0
(0)
u1
(0)
u2
(0)
u0
1.9302 0.05512 0.0245301
1.92314 -0.034491 0.0308497
2 0 2.85924e-032
1.9302 0.05512 0.0245301
x0
(0)
x1
(0)
x2
(0)
x0
(空白)
(空白)
(1)
x0
(1)
x1
(1)
x2
(1)
x0
y0
(1)
y1
(1)
y2
(1)
y0
(1)
u0
(1)
u1
(1)
u2
(1)
u0
......
......
···
···
···
···
···
···
1.8607 0.207852 0.0287241
1.94431 0.234342 2.64236e-032
1.90144 0.310047 2.69221e-032
1.8607 0.207852 0.0287241
x0 E
(N −1)
x1 E
(N −1)
x2 E
(N −1)
x0 E
1.90144
1.99378
1.97516
1.90144
(0)
0.132877 0.0269012
0.0788166 2.7483e-032
0.157103 2.64662e-032
0.132877 0.0269012
(N −1)
(0)
(1)
(N −1)
y0 E
(N −1)
y1 E
(N −1)
y2 E
(N −1)
y0 E
(N −1)
u0 E
(N −1)
u1 E
(N −1)
u2 E
(N −1)
u0 E
図 3.1 （左）example2.plt, （右）example2.plt の意味．ただし，2.7483e-032=
2.7483 × 10−32 など.
x0 y0 z0
x1 y1 z1
P0
x2 y2 z2
x3 y3 z3
x4 y4 z4
P4
P2
P1
P3
図 3.2 （左）データファイル（右）描画結果
もう少し見栄えを良くするために，gnuplot のカラーマップを利用する．この場合は，
gnuplot2 のように入力すると，図 3.3（右）が出力される．なお，set で記述されるコマ
ンドは，(一回の起動中に) 一度入力すればよい．gnuplot2 では，カラーマップの利用を
宣言し（1 行目），色の勾配を指定している（2 行目）．また，xy 平面を図の底面に合わ
せている（3 行目）．
gnuplot 2:
12
gnuplot> set pm3d
gnuplot> set palette rgbformulae 33,13,10
gnuplot> set ticslevel 0
gnuplot> splot "example2.plt" with lines pal
"example2.plt"
"example2.plt"
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0.4
0.35
0.3
0.25
0.2
0.15
0.1
1
-1.5
-1
0
-2
0
-0.5
0
0.5
-0.5
1
1.5
1
0.05
0.5
-2
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0.45
0.5
-1.5
-1
0
-0.5
0
0.5
2 -1
-0.5
1
1.5
2 -1
図 3.3 （左）gnuplot1, （右）gnuplot2
最後に，表示した絵を EPS 形式のファイルに保存するには，gnuplot 3 のようにする．
set term は，出力先の指定，set output は出力ファイルの指定である．replot で再
描画をしているが，splot ... をもう一度書いても同じである．最後に，画の出力先を，
このコンピュータの画面に戻している．また，結果を PDF 形式のファイルで保存する
場合は，gnuplot 4 のようにする．
gnuplot 3:
gnuplot> set term postscript enhanced color
gnuplot> set output "example2.eps"
gnuplot> replot
gnuplot> set term win
(Windows の場合)
gnuplot> set term x11
(Linux, Mac (X Window) の場合)
gnuplot> set term aqua
(Mac (AquaTerm) の場合)
gnuplot 4:
gnuplot> set term pdf
gnuplot> set output "example2.pdf"
gnuplot> replot
問題 3.1. example2.edp を作成し，実行せよ．example2.plt ができていることを確認
して，gnuplot で可視化せよ．
問題 3.2. 色の指定の際に，次のようにするとどうなるかを試してみよ．
set palette rgbformulae 7,5,15
13
set palette rgbformulae 22,13,-31
set palette rgbformulae 23,28,3
set palette rgbformulae 3,11,6
注意 3.1. 計算終了時に必ず「停止」を押すこと．そうしないと，データの出力が完了し
ないことがある．
14
4 領域と三角形分割
この節の例では，偏微分方程式の計算はしない．
例 4.1 (example3.edp). 図 4.1（左）で示されている領域をパラメータ表示し．三角形
分割を作る．
例えば，次のようにする:
Γ1 :
Γ2 :
Γ3 :
Γ4 :
( )
( )
( )
x
4
4
=t
+ (1 − t)
y
4
−4
( )
( )
( )
x
−4
4
=t
+ (1 − t)
y
4
4
( )
( )
( )
x
−4
−4
=t
+ (1 − t)
y
−4
4
( )
( )
( )
x
4
−4
=t
+ (1 − t)
y
−4
−4
(t : 0 → 1),
(t : 0 → 1),
(t : 0 → 1),
(t : 0 → 1)
パラメータの正の向きを，左回りにとっていることに注意せよ．プログラム exam-
ple3.edp では，各辺上に (等間隔に)n 個の節点を取るようにしている．Freefem++
は，境界上の節点の配置に合わせて要素を作成する．結果は，図 4.1（右）の通り．
buildsmesh では，Γ1 → Γ2 → Γ3 → Γ4 のように，閉曲線を描く順に定義しているこ
とに注意せよ．Γ1 → Γ3 → Γ2 → Γ4 では曲線が閉じないので許容されない．
リスト 4.1
example3.edp
//example3.edp
int n = 20;
border G1(t = 0, 1) { x =
border G2(t = 0, 1) { x =
border G3(t = 0, 1) { x =
border G4(t = 0, 1) { x =
mesh Th = buildmesh(G1(n)
4; y = -4*(1-t)+4*t;}
4*(1-t)-4*t; y = 4;}
-4; y = 4*(1-t)-4*t;}
-4*(1-t)+4*t; y = -4;}
+ G2(n) + G3(n) + G4(n));
plot(Th, ps="example3.eps");
例 4.2 (example4.edp). 図 4.2 で示されている領域をパラメータ表示し．三角形分割を
作る．大きな領域から，小さな領域を「くり抜く」には，小さな領域の境界のパラメータ
表示を右回りにする．すなわち，Γ5 , . . . , Γ9 のパラメータ表示は，次のようになる（いず
れも，パラメータの正の向きが，右回りになっていることに注意せよ）．結果は，図 4.3
に示されている．
15
Γ2
(−4, 4)
(4, 4)
Γ1
Γ3
(−4,−4)
(4, −4)
Γ4
図 4.1
Γ2
(−4, 4)
(4, 4)
(−2, 3)
Γ9
Γ6
(−3, 2)
Γ8
(−1, 2)
(−2, 1)
Γ9
Γ1
Γ7
(−2, 3)
Γ6
(−3, 2)
Γ3
Γ5
Γ5
(−1, 2)
Γ8
(2, −2)
1.5
Γ7
(−2, 1)
(2, −2)
1.5
(−4,−4)
Γ4
(4, −4)
図 4.2
Γ5 :
Γ6 :
Γ7 :
Γ8 :
Γ9 :
( ) (
)
x
2 + 1.5 cos t
=
y
−2 − 1.5 sin(t)
( )
( )
( )
x
−1
−2
=t
+ (1 − t)
y
2
3
( )
( )
( )
x
−2
−1
=t
+ (1 − t)
y
1
2
( )
( )
( )
x
−3
−2
=t
+ (1 − t)
y
2
1
( )
( )
( )
x
−2
−3
=t
+ (1 − t)
y
3
2
16
(t : 0 → 2π),
(t : 0 → 1),
(t : 0 → 1),
(t : 0 → 1),
(t : 0 → 1).
リスト 4.2
example4.edp
//example4.edp
int n = 8;
border G1(t = 0, 1) { x = 4; y = -4*(1-t)+4*t;}
border G2(t = 0, 1) { x = 4*(1-t)-4*t; y = 4;}
border G3(t = 0, 1) { x = -4; y = 4*(1-t)-4*t;}
border G4(t = 0, 1) { x = -4*(1-t)+4*t; y = -4;}
border G5(t = 0, 2*pi) { x = 2 + 1.5*cos(t); y = -2 - 1.5*sin(t);}
border G6(t = 0, 1) { x = -t - (1-t)*2; y = 2*t + 3*(1 - t);}
border G7(t = 0, 1) { x = -2*t - (1-t); y = t + 2*(1 - t);}
border G8(t = 0, 1) { x = -3*t - (1-t)*2; y = 2*t + (1 - t);}
border G9(t = 0, 1) { x = -2*t - (1-t)*3; y = 3*t + 2*(1 - t);}
mesh Th = buildmesh(G1(4*n) + G2(4*n) + G3(4*n) + G4(4*n)
+ G5(6*n) + G6(n) + G7(n) + G8(n) + G9(n));
plot(Th, ps="example4.eps");
図 4.3
問題 4.1. 図 4.4 で示した領域を，パラメータ表示し，三角形分割を作れ．
図 4.4
17
5 例題集
例 5.1 (example5.edp). Ω = (0, 1) × (0, 1) で，
−2∆u + 3u = ex+y in Ω,
u = sin(2πx) on Γ1 ,
u = 0 on Γ2
を解く．ただし，Γ1 = {(x, 0) | 0 < x < 1}, Γ2 = ∂Ω\Γ とする．f (x, y) = ex+y ,
g(x, y) = sin(2πx) とおくと，弱形式は，
u ∈ H 1 (Ω),
)
∫∫
∫∫ (
∫∫
∂u ∂v
∂u ∂v
uv dxdy =
f v dxdy
2
+
dxdy + 3
∂x ∂x ∂y ∂y
Ω
Ω
Ω
u = g on Γ1 ,
u = 0 on Γ2
(∀v ∈ H01 (Ω)),
となる．
リスト 5.1
example5.edp
// example5.edp
func f = exp(x + y);
func g = sin(2*pi*x);
func g0 = 0;
int n = 60;
border G1(t = 0, 1) {x = t; y = 0;}
border G2(t = 0, 1) {x = 1; y = t;}
border G3(t = 0, 1) {x = 1 - t; y = 1;}
border G4(t = 0, 1) {x = 0; y = 1 - t;}
mesh Th = buildmesh(G1(n) + G2(n) + G3(n) + G4(n));
fespace Vh(Th, P1);
Vh u,v;
solve elliptic(u, v) =
int2d(Th)(2*dx(u)*dx(v) + 2*dy(u)*dy(v)) +
int2d(Th)(3*u*v) - int2d(Th)(f*v) + on(G1, u = g) + on(G2,G3,G4, u = g0);
plot(u, ps="example5.eps");
例 5.2 (example6.edp). Ω = (0, 1) × (0, 1) とし，Neumann 境界値問題
−2∆u + 3u = ex+y in Ω,
∂u
= 0 on ∂Ω
∂n
を解く．弱形式は，f (x, y) = ex+y とおいて，
u ∈ H 1 (Ω),
)
∫∫ (
∫∫
∫∫
∂u ∂v
∂u ∂v
2
+
dxdy + 3
uv dxdy =
f v dxdy
∂x ∂x ∂y ∂y
Ω
Ω
Ω
18
(∀v ∈ H 1 (Ω)).
リスト 5.2
example6.edp
// example6.edp
func f = exp(x + y);
int n = 60;
border G1(t = 0, 1) {x = t; y = 0;}
border G2(t = 0, 1) {x = 1; y = t;}
border G3(t = 0, 1) {x = 1 - t; y = 1;}
border G4(t = 0, 1) {x = 0; y = 1 - t;}
mesh Th = buildmesh(G1(n) + G2(n) + G3(n) + G4(n));
fespace Vh(Th, P1);
Vh u,v;
solve elliptic(u, v) =
int2d(Th)(2*dx(u)*dx(v) + 2*dy(u)*dy(v)) + int2d(Th)(3*u*v) - int2d(Th)(f*v);
plot(u, ps="example6.eps");
例 5.3 (example7.edp). Ω = {(x, y) | x2 +y 2 < 1}, Γ1 = {(x, y) | x2 +y 2 = 1, y > 0},
Γ2 = {(x, y) | x2 + y 2 = 1, y < 0} とおいて，混合境界値問題
−∆u = 1 in Ω,
u = 1 on Γ1 ,
∂u
= 0 on Γ2
∂n
を解く．弱形式は，
V = {v ∈ H 1 (Ω) | v = 0 on Γ1 },
f (x, y) = 1
とおいて，
∫∫ (
Ω
∂u ∂v
∂u ∂v
+
∂x ∂x ∂y ∂y
)
u∈V :
∫∫
dxdy =
f v dx
(∀v ∈ V ),
Ω
リスト 5.3
example7.edp
// example7.edp
func f = 1;
func g = 1;
int n = 30;
border G1(t = 0, pi) {x = cos(t); y = sin(t);}
border G2(t = 0, pi) {x = cos(t + pi); y = sin(t + pi);}
mesh Th = buildmesh(G1(n) + G2(n));
19
u = 1 on Γ1 .
fespace Vh(Th, P1);
Vh u,v;
solve elliptic(u, v) =
int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) - int2d(Th)(f*v) + on (G1, u = g);
plot(u, ps="example7.eps");
例 5.4 (example8.edp). 例 5.3 と同じ Ω, Γ1 , Γ2 で，
−∆u + x2 u = ex+y in Ω,
u = 0 on Γ1 ,
∂u
= 1 on Γ2
∂n
を解く．b(x, y) = x2 , f (x, y) = ex+y , g(x, y) = 1 とおくと，Green の公式より，
φ ∈ C ∞ (Ω) に対して，
)
∫∫
∫
∫∫
∫∫ (
∂u ∂v
∂u
∂u ∂v
buv dxdy =
f v dxdy
+
dxdy +
v dS +
∂x ∂x ∂y ∂y
Ω
∂Ω ∂n
Ω
Ω
なので，弱形式は，V = {v ∈ H 1 (Ω) | v = 0 on Γ1 } とおいて，
∫∫ (
Ω
∂u ∂v
∂u ∂v
+
∂x ∂x ∂y ∂y
)
∫∫
dxdy +
u ∈ V,
∫
∫∫
buv dxdy =
Ω
gv dS +
Γ2
Freefem++ では，
∫
int1d(Th, G2)(g2*v) ←
gv dS
Γ2
のように記述する．
リスト 5.4
example8.edp
// example8.edp
func f = exp(x + y);
func b = x*x;
func g1 = 0;
func g2 = 1;
int n = 30;
border G1(t = 0, pi) {x = cos(t); y = sin(t);}
border G2(t = 0, pi) {x = cos(t + pi); y = sin(t + pi);}
mesh Th = buildmesh(G1(n) + G2(n));
fespace Vh(Th, P1);
Vh u,v;
solve elliptic(u, v) =
int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) + int2d(Th)((b*u)*v)
20
f v dxdy
Ω
(∀v ∈ V ).
- int1d(Th, G2)(g2*v) - int2d(Th)(f*v) + on (G1, u = g1);
plot(u, ps="example8.eps");
問題 5.1. 上記の例題を作成し，実行せよ．
問題 5.2. 実行結果をファイルに書き出し，gnuplot で可視化せよ．
問題 5.3. 問題 4.1 で作成した Ω を用いて計算するように，例題のプログラムを変更し
て，実行せよ．Γ1 , Γ2 などは，適当に決めよ．
21
6 有限要素解の収束
Poisson 方程式の境界値問題
−∆u = f
in Ω,
u=0
on ∂Ω
の有限要素近似 uh の収束について復習しておこう．Ω を多角形領域，Th をその三角形
分割とし，Vh ⊂ H01 (Ω) を連続区分一次要素とする．このとき，Poisson 方程式の有限要
素近似は，
∫∫
∫∫
uh ∈ Vh ,
∇uh · ∇vh dxdy =
Ω
f vh dxdy
(∀vh ∈ Vh )
Ω
で与えれる．
{Th }h>0 が正則な三角形分割の族で，u ∈ H 2 (Ω) ならば，
∥∇u − ∇uh ∥L2 (Ω) ≤ Ch|u|H 2 (Ω)
を満たすような正定数 C が存在する．さらに，Ω が凸多角形領域ならば，
∥u − uh ∥L2 (Ω) ≤ C ′ h2 |u|H 2 (Ω)
を満たすような正定数 C ′ が存在する．
例として，
−∆u = 2π 2 sin(πx) sin(πy)
in Ω,
u=0
on ∂Ω
を Ω = (0, 1) × (0, 1) において考える．この問題の厳密解は，
u(x, y) = sin(πx) sin(πy)
となる．
次の量を観察しよう：
eh = ∥∇u − ∇uh ∥L2 (Ω) ,
ρh =
log e2h − log eh
,
log(2h) − log(h)
Eh = ∥u − uh ∥L2 (Ω) ,
Rh =
log E2h − log Eh
.
log(2h) − log(h)
図 6.1 に示す規則メッシュを用いる．error1.edp（リスト 6.1）にこの計算を行うプログ
ラムを示す*2 ．
結果を，図 6.2 と表 6.1 に示す．これらの結果から，
eh ≈ C1 h,
Eh ≈ C2 h2
(C1 , C2 適当な定数),
が観察でき，これらは理論的な結果と整合している．
*2
オリジナル版の error1.edp にあった間違いを，木村正人先生（金沢大学）に指摘して修正して頂きまし
た．[2015 年 5 月 8 日]
22
図 6.1
1
log ERROR
0.1
0.01
0.001
H1 error
L2 error
0.0001
0.01
0.1
1
mesh size h
図 6.2
誤差の挙動
h
eh
ρh
Eh
Rh
0.353553
0.838452
—
0.079064
—
0.176777
0.431591
0.96
0.021119
1.90
0.088388
0.217113
0.99
0.005364
1.98
0.044194
0.108122
1.01
0.001337
2.00
0.022097
0.052783
1.03
0.000325
2.04
表 6.1
eh , ρh , eh ，Rh
23
リスト 6.1 error1.edp
// error1.edp
func exact = sin(pi*x)*sin(pi*y); // exact solution
func f = 2.0*pi*pi*exact; // rifht-hand side function
func g = exact; // Dirichlet boundary condition
real hsize, hold; // mesh size
real errh1, errh1old, errl2, errl2old, rateh1, ratel2;
int n, nn;
// fine triangulation
nn = 256;
mesh Th2 = square(nn, nn);
fespace Vh2(Th2, P1);
Vh2 uproj, w, uex;
uex = exact;
// output file
ofstream f1("error1.dat");
// n=4,8,16,32,64
n = 2;
errh1old = errl2old = 1.0;
hold = 1.0;
for (int i = 1; i < 6; i++) {
n = 2*n;
mesh Th = square(n, n);
plot(Th, ps="error1.eps");
fespace Vh(Th, P1);
Vh u, v, hh = hTriangle;
hsize = hh[].max;
solve Poisson(u,v) =int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v)) - int2d(Th) ( f
*v ) + on(1,2,3,4,u = g);
// computation of error using the fine triangulation
uproj = u; // projection of u into the fine triangulation
w = uproj - uex; // error function
errh1 = sqrt( int2d(Th2)(dx(w)*dx(w) + dy(w)*dy(w)) ); // H1-error
errl2 = sqrt( int2d(Th2)(w^2) ); // L2-error
// computation of rates
rateh1 = (log(errh1) - log(errh1old))/(log(hsize) - log(hold));
ratel2 = (log(errl2) - log(errl2old))/(log(hsize) - log(hold));
errh1old = errh1;
errl2old = errl2;
hold = hsize;
// output results
f1 << hsize << " "
<< errh1 << " " << rateh1 << " "
<< errl2 << " " << ratel2<< " " << endl;
}
24
7 非線形問題への応用：Nehari の変分原理
p > 1 に対して，非線形楕円型方程式
− ∆u = up , u ≥ 0 in Ω,
u = 0 on Γ ≡ ∂Ω
(7.1)
を考えよう．Ω は R2 の有界領域である．Ω の形状によっては非自明な解は複数個有り
得る．それらは V = H01 (Ω) 上の汎関数
∫∫ (
J(v) =
Ω
1
1
2
p+1
|∇v| −
|v|
2
p+1
)
dxdy
の臨界点として特徴付けられる．その中でも，特に，
u∈V :
J(u) = inf J(v)
v∈N
を満たす u を Nehari の解と呼ぶことにしよう．ここで，
{
N = v∈V
∫∫
}
∫∫
2
p+1
|∇v| dxdy =
|v| dxdy
Ω
Ω
であり，これをを Nehari の多様体と呼ぶ．上記の変分構造に基づいて，Nehari の解を
求めるために，次の反復解法を考えることができる (cf. [2] の §6.1 や [6])．
• 初期関数 u0 (x) ≥ 0 を選ぶ．
∞
∞
• {wn (x)}∞
n=0 ⊂ V , {un (x)}n=1 ⊂ N , {λn }n=0 ⊂ R+ を次で求める：
)
∫∫ (
∫∫
∂wn ∂v
∂wn ∂v
wn ∈ V,
+
dxdy =
upn v dxdy, ∀v ∈ V ;
∂x ∂x
∂y ∂y
Ω
Ωh
(
)1/(p−1)
∫∫
∫∫
(1)
λn
(1)
2
(2)
|wn |p+1 dxdy;
λn =
λn =
|∇wn | dxdy, λn =
(2)
Ω
Ω
λn ,
un+1 = λn wn .
定理 7.1 ([2], §6.1).
• λn < 1 であり，単調に λn ↑ 1.
• {∥∇un ∥L2 (Ω) } は単調減少かつ (正の定数により) 下から有界．
• un+1 = T un で作用素 T : N → N を定めると，
J(T v) ≤ J(v)
(v ∈ N ).
さらに，v ∈ N が (7.1) の解であるときのみ等号が成立する．
このアルゴリズムを Freefem++ で実行するためのプログラムが nehari1.edp である．
25
リスト 7.1
nehari1.edp
//nehari1.edp
func g = 0;
func initial = 1;
real p = 3.0, q = 1.0/(p - 1);
real s1, s2, lambda, ep;
int i, j, iter, maxiter = 200, k = 20;
// output files
ofstream f1("nehari1.dat");
ofstream f2("nehari1.plt");
// domain
border G1(t = 0, 3) {x = t; y = 0;}
border G2(t = 0, pi/2) {x = 3*cos(t); y = 3*sin(t);}
border G3(t = 0, 3) {x = 0; y = 3 - t;}
border G4(t = 0, 2*pi) {x = 1.9 + 0.8*cos(t); y = 0.9 - 0.8*sin(t);}
border G5(t = 0, 2*pi) {x = 0.7 + 0.5*cos(t); y = 2.3 - 0.5*sin(t);}
// triangulation
mesh Th = buildmesh(G1(k)+G2(2*k)+G3(k)+G4(2*k)+G5(k));
// P1 element
fespace Vh(Th,P1);
Vh u,v,w,f, hh = hTriangle;
// n = 0
u = initial;
// n = 1, 2, ...
lambda = 100;
iter = 0;
ep = hh[].max*0.1;
while(abs(lambda - 1.0) > ep){
iter++;
if(iter > maxiter) break;
f = u^p;
solve poisson(w, v)=int2d(Th)(dx(w)*dx(v) + dy(w)*dy(v))
- int2d(Th)(f*v) + on(G1,G2,G3,G4,G5,w = g);
s1 = int2d(Th)(w^(p+1));
s2 = int2d(Th)(dx(w)*dx(w) + dy(w)*dy(w));
lambda = (s2/s1)^q;
u = lambda*w;
f1 << iter << " "<< lambda << " " << endl;
}
//results
f2 << "# iteration:" << iter << endl << endl;
for(i = 0; i < Th.nt ; i++){
for(j = 0; j < 3; j++){
f2 << Th[i][j].x << " "<< Th[i][j].y << " " << u[][Vh(i,j)] << endl;
}
26
f2 << Th[i][0].x << " " << Th[i][0].y << " " << u[][Vh(i,0)] << endl <<
endl << endl;
}
//
plot(Th, ps="nehari1mesh.eps");
nehari1.edp を詳しく見ていこう．
//nehari1.edp
func g = 0;
func initial = 1;
real p = 3.0, q = 1.0/(p - 1);
real s1, s2, lambda, ep;
int i, j, iter, maxiter = 200, k = 60;
ここでは，使用する変数を定義している．real で定義したものは，（倍精度）実数型と
なる．ep は反復の停止条件，maxiter は最大反復回数である．
// output files
ofstream f1("nehari1.dat");
ofstream f2("nehari1.plt");
// domain
.....
// triangulation
mesh Th = buildmesh(G1(k)+G2(2*k)+G3(k)+G4(2*k)+G5(k));
データ出力用に，2 つのファイル nehari1.dat と nehari1.plt を用意する．
nehari1.plt は，example2.edp と同様に，gnuplot 出力用のデータを保存する．一
方で，nehari1.dat には，各反復における λn の値を出力しておく．これで，λn → 1
となっているかどうかを確認できる．
// P1 element
fespace Vh(Th,P1);
Vh u,v,w,f, hh = hTriangle;
27
連続区分一次要素を定義する．hh=hTriangle は，各要素の最大辺長を表す関数である．
// n = 0
u = initial;
u0 (x) ≡ 1 とする．
// n = 1, 2, ...
lambda = 100;
iter = 0;
ep = hh[].max*0.1;
while(abs(lambda - 1.0) > ep){
iter++;
if(iter > maxiter) break;
f = u^p;
solve poisson(w, v)= ....
s1 = int2d(Th)(w^(p+1));
s2 = int2d(Th)(dx(w)*dx(w) + dy(w)*dy(w));
lambda = (s2/s1)^q;
u = lambda*w;
f1 << iter << " "<< lambda << " " << endl;
}
while(abs(lambda - 1.0) > ep){ ... }は，abs(lambda-1.0)>ep が正しいとき
にのみ，{ ... }を実行し，そうでなければ実行せずに先に進む，という命令である．
abs(x) は x の絶対値．いまは，|λn − 1| = abs(lambda-1.0) ≤ ε となったら，収束し
たとみなし，計算を終了するようになっている．ただし，ε = ep は，メッシュの細かさ
（≈ 近似解の精度）に合わせて，ep = hh[].max*0.1 と定義している．hh[].max は，
Th = Th の最大辺長を表す．ただし，これだけでは，反復列が収束しなかったときに，
永久に計算が実行されることになるので，最大反復回数 maxiter を定義しておき，反復
回数 iter が maxiter を超えたら，強制的に反復計算をやめる様にしてある．それの命
令が，if(iter > maxiter) break; である．なお，iter は（現在の）反復回数を表わ
し，iter++ は iter=iter+1 の意味である．
28
//results
f2 << "# iteration:" << iter << endl << endl;
for(i = 0; i < Th.nt ; i++){
....
}
この部分は，example2.edp と同じ．
例 7.1. nehari1.edp を k=20, 40, 60, 80 に対して実行する．それぞれの実行結果を図
7.1—7.3 に示した．また，図 7.4 に，λn の収束履歴を示した．
6
5
6
4
5
3
4
2
1
3
0
2
3
1
0
0
2.5
2
0.5
1.5
1
1.5
1
2
2.5
0.5
3 0
図 7.1 nehari1.edp (k = 20)
5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0
5
4.5
4
3.5
3
2.5
2
1.5
3
1
2.5
0.5
00
2
0.5
1.5
1
1.5
1
2
2.5
0.5
3 0
図 7.2 nehari1.edp (k = 40)
問題 7.1. nehari1.edp を作成・実行し，結果を gnuplot で可視化せよ．メッシュの細か
さや，初期関数 initial をいろいろと変えてみよ．計算の都度，nehari1.dat を見て，
収束を確認すること．
問題 7.2. 問題 4.1 で作成した領域等を用いて，問題 7.1 と同様の実験・考察を行なえ．
29
5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0
5
4.5
4
3.5
3
2.5
2
1.5
3
1
2.5
0.5
2
00
1.5
0.5
1
1.5
1
2
2.5
0.5
3 0
図 7.3 nehari1.edp (k = 60)
1
0.95
0.9
lambda_n
0.85
0.8
0.75
0.7
0.65
0.6
k=20
k=40
k=60
0.55
2
2.5
3
3.5
4
4.5
5
iteration
図 7.4
横軸: 反復回数 n，縦軸: λn (nehari1.edp, n = 20, 40, 60, 80)
30
8 非線形問題への応用：スケーリング反復法
p > 1 に対して，非線形楕円型方程式
− ∆u + au = bup , u ≥ 0 in Ω,
u = 0 on Γ ≡ ∂Ω
(8.1)
を考えよう．Ω は R2 の有界領域，a ≥ 0, b > 0 は定数である．
v ∈ V = H01 (Ω) と α > 0 が，
− ∆v + av = bαv p , v ≥ 0 in Ω,
v = 0 on Γ,
∥v∥∞ = ∥v∥L∞ (Ω) = 1
(8.2)
を満たしているとしよう．このとき，β = α1/(p−1) > 0 と定義して，u = βv とおくと，こ
の u は (8.1) を満たしている．スケーリング反復法 (SIA, Scaling iterative algorithm)
は，(8.2) の解 v を反復的に求めることで，(8.1) の解 u の近似値を得る方法である．
■ スケーリング反復法 (SIA, Scaling iterative algorithm; [1] の §2.2):
1
u0 (x) とする.
β0
∞
∞
∞
∞
• {wn (x)}∞
n=1 , {vn (x)}n=1 , {un (x)}n=1 , {αn }n=1 , {βn }n=1 を次で生成:
• 初期値 u0 (x) を選び，β0 = ∥u0 ∥∞ , v0 (x) =
p
, wn ≥ 0 in Ω, wn = 0 on Γ;
−∆wn + awn = bvn−1
1
αn =
;
vn (x) = αn wn (x);
∥wn ∥∞
βn = αn1/(p−1) ;
un (x) = βn vn (x).
これを実行するプログラムが，sia1.edp である．なお，反復の停止条件は，相対的な増分
δn =
∥∇vn − ∇vn−1 ∥
∥∇vn−1 ∥
をチェックし，δn < ε となった時点で収束とみなし，un = βn vn = beta*v を出力して
いる．また，初期値 u0 (x, y) は，u0 (x, y) = x2 + y 2 としている．
リスト 8.1
sia1.edp
//sia1.edp
func g = 0;
func initial = x*x + y*y;
real a = 0.0, b = 1.0, p = 3.0;
real beta, alpha, delta, norm, ep = 1.0e-4;
int i, j, iter, maxiter = 100, k = 20;
// output files
ofstream f1("sia1.dat");
ofstream f2("sia1.plt");
31
// domain
border G1(t = 0, 3) {x = t; y = 0;}
border G2(t = 0, pi/2) {x = 3*cos(t); y = 3*sin(t);}
border G3(t = 0, 3) {x = 0; y = 3 - t;}
border G4(t = 0, 2*pi) {x = 1.9 + 0.8*cos(t); y = 0.9 - 0.8*sin(t);}
border G5(t = 0, 2*pi) {x = 0.7 + 0.5*cos(t); y = 2.3 - 0.5*sin(t);}
// triangulation
mesh Th = buildmesh(G1(k)+G2(2*k)+G3(k)+G4(2*k)+G5(k));
// P1 element
fespace Vh(Th,P1);
Vh u,v,vv,w,f;
// n = 0
u = initial;
beta = u[].max;
vv = v = u/beta;
// n = 1, 2, ...
delta = 100;
iter = 0;
while(delta > ep){
iter++;
if(iter > maxiter) break;
f = b*v^p;
solve poisson(u, v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v) + a*u*v)
- int2d(Th)(f*v) + on(G1,G2,G3,G4,G5,u = g);
alpha = 1/(u[].max);
v = alpha*u;
w = v - vv;
norm = sqrt(int2d(Th)(dx(v)*dx(v) + dy(v)*dy(v)));
delta = sqrt(int2d(Th)(dx(w)*dx(w) + dy(w)*dy(w)))/norm;
vv = v;
f1 << iter << " "<< delta << " " << norm << endl;
}
//
beta = alpha^(1/(p-1));
u = beta*v;
//results
f2 << "# iteration:" << iter << endl << endl;
for(i = 0; i < Th.nt ; i++){
for(j = 0; j < 3; j++){
f2 << Th[i][j].x << " "<< Th[i][j].y << " " << u[][Vh(i,j)] << endl;
}
f2 << Th[i][0].x << " " << Th[i][0].y << " " << u[][Vh(i,0)] << endl <<
endl << endl;
}
//
plot(Th,wait=true, ps="sia1mesh.eps");
32
8
7
6
5
4
3
2
1
0
8
7
6
5
4
3
2
3
1
00
図 8.1
2.5
2
0.5
1.5
1
1.5
1
2
2.5
0.5
3 0
sia1.edp (k = 20)
8
7
6
5
4
3
2
1
0
8
7
6
5
4
3
2
3
1
00
図 8.2
2.5
2
0.5
1.5
1
1.5
1
2
2.5
0.5
3 0
sia1.edp (k = 40)
8
7
6
5
4
3
2
1
0
8
7
6
5
4
3
2
3
1
00
図 8.3
2.5
2
0.5
1.5
1
1.5
1
2
2.5
0.5
3 0
sia1.edp (k = 60)
例 8.1. sia1.edp を k = 20, 60 に対して実行した結果が，図 8.1, 8.2 である．一方，図
8.4 には，δn の収束の履歴を両対数目盛で表示した．
問題 8.1. sia1.edp を作成・実行し，結果を gnuplot で可視化せよ．メッシュの細かさ
や，初期関数をいろいろと変えてみよ．計算の都度，sia1.dat を見て，収束を確認する
こと．
問題 8.2. 問題 4.1 で作成した領域等を用いて，問題 8.1 と同様の実験・考察を行なえ．
問題 8.3. 例 7.1 と例 8.1 は同じ方程式を計算しているが，全く異なった解が得られてい
33
10
increment: log(delta_n)
1
0.1
0.01
0.001
0.0001
k=20
k=40
k=60
1e-05
1
10
100
iteration: log(n)
図 8.4
横軸: 反復回数 log n，縦軸: log δn (parabolic1.edp, k = 20, 60)
る．何が異なっているのだろうか？同じ解を捉えることはできるだろうか？
34
35
参考文献
[1] G. Chen, J. Zhou and W. M. Ni: Algorithms and visualization for solutions
of nonlinear elliptic equations, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 10
(2000) 1565–1612.
[2] H. Fujita, N. Saito and T. Suzuki: Operator Theory and Numerical Methods,
North-Holland, Elsevier, Amsterdam, 2001.
[3] F. Hecht, O. Pironneau, A. Le Hyaric and K. Ohtsuka: Freefem++ (Ver. 3.37),
http://www.freefem.org/ff++/, Laboratoire Jacques-Louis Lions (LJLL),
University of Paris VI.
[4] 大塚厚二，Freefem++ による有限要素解析入門，http://comfos.org/jp/
ffempp/index.html
[5] 大塚厚二，高石武史：有限要素法で学ぶ現象と数理—Freefem++ 数理思考プログ
ラミング—，共立出版，2014 年
[6] 鈴木貴，上岡友紀：偏微分方程式講義—半線形楕円型方程式入門，培風館，2005 年．
36

Download Report