FreeFem の紹介
桂田 祐史

2007 年 9 月 26 日〜, 2016 年 2 月 14 日, 2021 年 11 月 7 日

目次
1

手短な紹介

1

2

百聞は一見に如かず
2.1 Mac に最新版をインストールする . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 簡単な例題で超特急の解説 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
3

3

おまけ: gnuplot で可視化

8

A 参考: FreeFem++ で菊地 [1] の Poisson 方程式の例題を解く

9

B 2 次元 Laplacian の固有値問題を FreeFem++ で解く

13

C FreeFem++ のマニュアル §3.9 から: Stokes 方程式の cavity ﬂow

18

D 熱方程式に対する後退 Euler 法
この文書の最新版は

21

http://nalab.mind.meiji.ac.jp/~mk/labo/text/welcome-to-freefem/
または

http://nalab.mind.meiji.ac.jp/~mk/labo/text/welcome-to-freefem.pdf
で閲覧可能です。
FreeFem の文法については、一新された FreeFem Documentation1 を読むのが本筋だと思います。(私は 「FreeFem++ノート」2 というメモも書いてありましたが、今後どうしようか
思案中です。)

手短な紹介

1

偏微分方程式の汎用数値解法には、古くから差分法 (ﬁnite diﬀerence method, FDM) と有限要素法 (ﬁnite element method, FEM) という二大手法があります (有限体積法 (ﬁnite volume
method) も良く使われるので、「二大」は止めた方が良いでしょうか…)。

1
2

https://doc.freefem.org/documentation/
http://nalab.mind.meiji.ac.jp/~mk/labo/text/freefem-note.pdf

1

有限要素法は、解析学の世界でスタンダードとなっている弱解の方法を数値計算に応用したもので、微分方程式の弱定式化の議論とほぼ並行した形で近似解の議論がなされ、数学者になじみやすいものです。また有限要素法は、プログラムの自動生成がしやすいという特徴を持っています3 。このため、従来から、プログラムを自動生成するシステムを開発する試みが色々なされてきました。

FreeFem4 は (FreeFem++5 というのは、不正確な言い方かもしれませんが、とりあえず “
一段前のバージョン” と呼んでおきます)、パリ第 6 大学 J. L. Lions 研究所の Frédéric Hecht,
Oliver Pironneau, A. Le Hyaric, 広島国際学院大学の大塚厚二氏らが開発した、2 次元, 3 次元問題を有限要素法で解くための、一種の PSE (problem solving environment) で、ソースコー
ドとマニュアル (幸い英文)、主なプラットホーム (Windows, Mac, Linux) 向けの実行形式パッ
ケージが公開されています。メインの開発者による解説 [?] を紹介しておきます。
マニュアル PDF (700 ページ超) が、FreeFEM Documentation6 という URL で入手できます (この URL は固定されたでしょうか…)
現在は WWW 上のドキュメンテーション FreeFem++ Documentation7 を参照するのが普通になるのかもしれません。
FreeFem では、ユーザーが境界曲線を指定することで定義した 2 次元領域に対して、自動的に三角形要素分割を生成し、弱形式により記述された問題に対して、種々の有限要素 (関数)空間を用いて弱解を求め、可視化することが可能になっています (最近では 3 次元領域の問題も解けるようになっているようですが、筆者はそれについてほとんど知らないので、言及するのを止めておきます。)。
問題の領域や弱形式などは専用の言語で記述した短いプログラム (ものによっては 20〜30
行程度) として与え、それを実行することでシミュレーションを行います。
すべての問題がこのソフトで扱えるわけではありませんが、扱える問題は結構幅広く、それが出来るときは多くの場合、C 言語のようなプログラミング言語で一からプログラムを書くのに比べて、桁違いに短い時間でシミュレーション実行までたどり着けます。(プログラムの行数自体は、2 桁とは行きませんが、1.5 桁くらい小さいように思います。)
大塚・高石「有限要素法で学ぶ現象と数理―FreeFem++」共立出版 (2014) という解説書が出版され、そのサポートサイト「有限要素法で学ぶ現象と数理」8 が出来ました。有益な情報が日本語で得られるようになりました。
日本応用数理学会で、FreeFem++ のセミナーが何回か開かれています。私が参加できた直近 2 回の情報は• ソフトウェアセミナー：FreeFem++による有限要素プログラミング −上級編−9 (2016/6/4,5)
資料 https://www.ljll.math.upmc.fr/~suzukia/FreeFempp-tutorial-JSIAM2016b/
index.html• ソフトウェアセミナー：FreeFem++で学ぶ現象と数理10 (2016/2/11, 12)
資料 http://www.ljll.math.upmc.fr/~suzukia/FreeFempp-tutorial-JSIAM2016
3
こう書くと差分法のシンパに反発されそうな気もしますが、野次馬の観察 (岡目八目) によると、有限要素法に一日の長があるのは否定できないと思います。
4https://freefem.org/
5
http://www3.freefem.org/
6
https://doc.freefem.org/pdf/FreeFEM-documentation.pdf
7
https://doc.freefem.org/introduction/index.html
8
http://comfos.org/jp/ffempp/book/9
http://na.cs.tsukuba.ac.jp/acmi/?p=403
10
http://na.cs.tsukuba.ac.jp/acmi/?p=383

2

幸い、2 回とも都合がついて参加でき、有益な情報が得られました。関係者のご尽力に感謝いたします。
講師の鈴木厚先生には、2007 年の九州大学で行われたチュートリアルでもお世話になりました。そのとき頂いた資料をネタにして、自分の学生に FreeFem++ を教えて来ましたが、今後は、上のセミナーの資料を咀嚼して、学生に教えていくことになると思います。

2

百聞は一見に如かず

2.1

Mac に最新版をインストールする

FreeFem++ 本体をインストールするには、WWW サイト Freefem++ Home Page11 から、FreeFem++-4.9-full-MacOS 10.11.pkg のようなインストーラー (4.9 はバージョン番号
で、これは当然変化します。また MacOS 10.11 は macOS 10.11 以降に対応する、という意
味です。) を入手して実行します (Ctrl+クリックして「開く」)。
例えば version 4.9 の場合は、主な実行可能プログラムは、/usr/local/ff++/mpich3/binというディレクトリィに置かれています (/Application ディレクトリィにも FreeFem++.app
が置かれますが、そのサイズは小さく、大部分は /usr/local/ff++/ の下にあります)。
ターミナルから実行するには、環境変数 PATH が設定されている必要がありますが、最近の FreeFem++ では、/etc/paths.d/FreeFem++ というファイルを用意することによって “自動的に” 設定されます (中身は /usr/local/ff++/mpich3/bin という 1 行です)。ですから特に意識する必要はないでしょう。
(新しいバージョンをインストールしてうまく行かない場合は、/usr/local/ff++ の下を探検してみることを勧めます12 。)「起動できない」現象数理学科の学生に: 時々シェルの設定ファイル (~/.profile とか~/.bash profileとか) がおかしくなっていて、/etc/paths.d/FreeFem++ での設定が有効にならないケースを見かけます (anaconda がらみの問題？)。そういう人は、設定ファイルの末尾に

export PATH=$PATH:/usr/local/ff++/mpich3/binという 1 行を書き加えてから、シェルを起動して (ターミナルで新しいウィンドウを出すとか、
ターミナルそのものを再起動する) みて下さい。
かつては、GUI でプログラムの開発＆実行の出来る FreeFem++-cs という IDE (統合開発環境) があって、以前は私も身の回りの学生に推奨していましたが、最近ではオリジナルのFreeFem++ だけで十分と考えるようになりました。
何か適当なテキスト・エディター (Atom (atom), Visual Code Studio (code), Emacs (emacs)など何でも) を選んで、テキスト・エディターでプログラムを入力・編集して、ターミナルでFreeFeem++ コマンドを使って実行する、というやり方を勧めます。

2.2

簡単な例題で超特急の解説

Poisson 方程式の境界値問題で説明します (有限要素法を自習したい場合には、菊地文雄,『有限要素法概説』, サイエンス社, という参考書がイチオシです, 桂田研の学生で希望する人は申し出て下さい)。
1112

http://www.freefem.org/
ターミナルで、open /usr/local/ff++ とすると、Finder で が開かれます。

3

R2 内の有界な領域 (連結開集合) Ω における Poisson 方程式の Dirichlet 問題(PE)− △u = f(DBC)u = 0 (on Γ := ∂Ω)(in Ω),

を考えます (弱解の方法の例題として、もっとも簡単な、定番の問題です)。
Poisson 方程式 (PE) に、任意の v ∈ C0∞ (Ω) (C0∞ (Ω) は、Ω 内に compact な台を持つ C ∞級の関数全体の集合) をかけて、Green の積分公式 (多次元版部分積分) を用いると、次式が得られます。
ZZ ZZ∂u ∂v ∂u ∂v(1)+
dx dy =f (x, y)v(x, y) dx dy (v ∈ C0∞ (Ω)).
∂x ∂x ∂y ∂yΩΩ
逆に u = 0 (on Γ) を満たす滑らかな u が、この条件 (1) を満たせば、u は Poisson 方程式(PE) を満たすことも容易に分かります。
微分方程式 (PE) を満たす u を求める代りに、(1) を満たす u を見出すことを考えましょう。
解の存在を容易に保証できるようにするため、関数空間を完備化するのが便利です (解を極限の論法で得よう、という解析学の常套手段を使うため)。具体的には、以下に定義する H01 (Ω)というソボレフ空間 (一般化された導関数を持つ関数からなる関数空間の一種) で u を探すことにします。
(PE),(DBC) の弱解の定義

u ∈ H01 (Ω) かつ
ZZ 
Ω



∂u ∂v ∂u ∂v+∂x ∂x ∂y ∂yZZdx dy =

H (Ω) =∂u∈ L2 (Ω) (j = 1, 2)
u ∈ L (Ω)
∂xj

ZZ

ZZ

1

(v ∈ H01 (Ω)).

Ω


H 1 (Ω) とは0

f (x, y)v(x, y) dx dy2という関数の集合に(u, v)H 1 :=∇u(x, y) · ∇v(x, y) dx dyu(x, y)v(x, y) dx dy +ΩΩという内積を導入すると、H 1 (Ω) は Hilbert 空間 (完備な内積空間) になる。

H01 (Ω) := C0∞ (Ω) の H 1 (Ω) での閉包とおく。H01 (Ω) の要素は、H 1 (Ω) の要素のうちで、ある意味で Γ 上 0 となるものとみなすことができる。
Ω を三角形の「整った」合併 Th で近似し、Th で連続で、各三角形上で 1 次関数であるような関数全体を Vh とします。Vh は H 1 (Ω) の有限次元部分空間とみなせ、グラフが角錐であるような関数 φ1 , . . . , φm を基底に取って、各元をその基底の線型結合として表現することが
できます:
)
( m
X
cj φj (cj ) ∈ Rm .
Vh =
j=1

4


(PE),(DBC) の弱解 u の定義 (書き換え)u ∈ H 1 (Ω) かつ

ZZ 
ZZ
∂u ∂v ∂u ∂v+
dx dy =f (x, y)v(x, y) dx dy∂x ∂x ∂y ∂yΩΩ(v ∈ H 1 (Ω), v = 0 on Γ),

u = 0 (on Γ).

有限要素解 uh の定義uh ∈ Vh かつ

ZZ 
ZZ
∂uh ∂vh ∂uh ∂vh+dx dy =f (x, y)vh (x, y) dx dy
∂x ∂x
∂y ∂y
Th
Th

(vh ∈ Vh , vh = 0 on Γ),

uh = 0 (on Γ).


次のようなプログラム poisson.edp を用意します (テキスト・エディターを使って作成します。現象数理学科学生は、mi, emacs, Atom などを習っているはずです。)。
poisson.edp13

// poisson.edp// http://nalab.mind.meiji.ac.jp/~mk/program/fem/poisson.edp// http://nalab.mind.meiji.ac.jp/~mk/labo/text/welcome-to-freefem/node4.html// 境界の定義 (単位円), いわゆる正の向き
border Gamma(t=0,2*pi) { x=cos(t); y=sin(t); }// 三角形要素分割を生成 (境界を 50 に分割)
mesh Th = buildmesh(Gamma(50));plot(Th,wait=true); // plot(Th,wait=true,ps="Th.eps");// 有限要素空間は P1 (区分的１次多項式) 要素
real [int] levels =-0.012:0.001:0.012;
fespace Vh(Th,P1);
Vh u,v;// Poisson 方程式 -△ u=f の右辺
func f = x*y;// 現在時刻をメモ
real start = clock();// 問題を解く
solve Poisson(u,v)= int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(f*v)+on(Gamma,u=0);// 可視化 (等高線)
plot(u,wait=true);//plot(u,viso=levels,fill=true,wait=true);// 可視化 (3 次元) --- マウスで使って動かせる//plot(u,dim=3,wait=true);plot(u,dim=3,viso=levels,fill=true,wait=true);// 計算時間を表示
cout << " CPU time= " << clock() - start << endl;ターミナルで poisson.edp を入手して実行 ($ はプロンプトなので入力しない)$ curl -O http://nalab.mind.meiji.ac.jp/~mk/program/fem/poisson.edp$ FreeFem++ poisson.edp


このプログラムを実行すると、円盤領域の三角形分割を表示したところで停止します (plot()に wait=true としてあるため)。先に進めるには enter (return) を入力します。

5

等高線が表示されたところでも、先に進めるには enter (return) を入力します。
最後にグラフを描画しますが、そのウィンドウを閉じるには esc を入力します。
キーボードからのコマンド[Enter]
[ESC]
+, Z
−, z
=
3
c
C
H, h
r
f
v
?



次のプロットへ
グラフィックスを閉じる
ズーム（イン）する
ズームアウトする
ズームの状態を元に戻す
2 次元/3 次元 表示の切り替え
ベクトルの矢印の大きさを短くする
ベクトルの矢印の大きさを長くする
z 方向のスケールを変える
リフレッシュする
カラーを塗りつぶすかどうか

help





(最近コマンドが増えています。ぜひ自分で ? として、色々試して下さい。)

図 1: v4.9 のコマンド一覧

plot(u,wait=true); のところを plot(u,wait=true,ps="poisson-disk.eps"); のようにすると、画面に表示するだけでなく、PostScript 形式でファイル出力できます (TEX 文書への取り込みが簡単です)。

6図 2: 円盤領域で Poisson 方程式 − △ u(x, y) = xy の同次 Dirichlet 問題を解く (三角形要素
分割と解の等高線)

図 3: dim=3 でグラフの鳥瞰図を描く

7

おまけ: gnuplot で可視化

3

従来、FreeFem++ では、等高線描画やベクトル場の表示などは出来るが、グラフの鳥瞰図
描画などはサポートしていなかった (今は plot(...,dim=3,...); とするだけで出来る)。そこで gnuplot を使ってグラフを描く方法を紹介する。このやり方は色々応用が効くので知っていて損はない。
poisson-g.edp14

// poisson-g.edp// 境界の定義 (単位円), いわゆる正の向き
border Gamma(t=0,2*pi) { x=cos(t); y=sin(t); }// 三角形要素分割を生成 (境界を 50 に分割)
mesh Th = buildmesh(Gamma(50));plot(Th,ps="Th.eps",wait=true);// 有限要素空間は P1 (区分的１次多項式) 要素
fespace Vh(Th,P1);
Vh u,v;// Poisson 方程式 -△ u=f の右辺
func f = x*y;// 現在時刻を記録
real start = clock();// 問題を解く
solve Poisson(u,v)= int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(f*v)+on(Gamma,u=0);// 可視化plot(u,ps="poisson-disk.eps");// 計算結果をファイルに出力{
ofstream ug("u-g.txt");
for (int i=0; i<Th.nt; i++) {for (int j=0; j<3; j++) {ug << Th[i][j].x << " " << Th[i][j].y << " " << u[][Vh(i,j)]<<endl;}
ug << Th[i][0].x << " " << Th[i][0].y << " " << u[][Vh(i,0)]<<"\n\n\n";}}// 計算時間を表示
cout << " CPU time= " << clock() - start << endl;


入手と実行




curl -O http://nalab.mind.meiji.ac.jp/~mk/program/fem/poisson-g.edp
FreeFem++ poisson-g.edp


実行が終わると “u-g.txt” というファイルが作成される。“u-g.txt” は次のようなデータ
である (3 次元空間における三角形が並んでいる)。

8
u-g.txt15

0.84068 0.327399 0.00414118
0.794983 0.411712 0.00532531
0.751583 0.286674 0.00623707
0.84068 0.327399 0.00414118

0.870933 0.153872 0.00249565
0.84068 0.327399 0.00414118
0.751583 0.286674 0.00623707
0.870933 0.153872 0.00249565(中略)......
0.1092 0.735914 0.00288613
0.140979 0.864349 0.00225845
0.0271079 0.868281 0.000433488
0.1092 0.735914 0.00288613-0.0867634 0.872212 -0.00146171-0.203275 0.85648 -0.00323151-0.110736 0.746398 -0.00290419-0.0867634 0.872212 -0.00146171


このデータを gnuplot で可視化するには、例えば次のようにする (図 4)。
poisson-g.g

set hidden3d
set palette rgbformulae 33,13,10
splot "u-g.txt" with lines palette

この内容を “gnuplot> ” という gnuplot のプロンプロトに入力すれば良いが、自動的にやるには、以下の 2 つのコマンドをターミナルで実行しても良い。
curl -O http://nalab.mind.meiji.ac.jp/~mk/program/fem/poisson-g.g
gnuplot poisson.g




参考文献
[1] 菊地文雄：有限要素法概説, サイエンス社 (1980), 新訂版 1999.

A

参考: FreeFem++ で菊地 [1] の Poisson 方程式の例題を解く

菊地 [1] に載っている Poisson 方程式の例題

(2)
(3)
(4)

−△u = fin Ω,

u = g1 in Γ1 ,∂u
= g2 in Γ2∂n

9

"u-g.txt"

0.015
0.01

0.015

0.005

0.01

0

0.005
0

-0.005-0.005

-0.01

-0.01

-0.015-0.015

-1
-0.5
0
0.5
1 -1

0
-0.2
-0.4
-0.6
-0.8

1
0.8
0.6
0.4
0.2

図 4: Poisson 方程式の解のグラフ表示 (gnuplot による)

(ただし、Ω = (0, 1) × (0, 1), Γ1 = {0} × [0, 1] ∪ [0, 1] × {0}, Γ2 = {1} × (0, 1] ∪ (0, 1] × {1},f = 1, g1 = 0, g2 = 0) を FreeFem++ を用いて解くにはどうするか。
例えば次のようなプログラムで解ける。

10
poisson-kikuchi.edp16

// poisson-kikuchi.edp// http://nalab.mind.meiji.ac.jp/~mk/program/fem/poisson-kikuchi.edp// 菊地文雄, 有限要素法概説, サイエンス社
int Gamma1=1, Gamma2=2;
border Gamma10(t=0,1) { x=0;
y=1-t; label=Gamma1; }border Gamma11(t=0,1) { x=t;
y=0;label=Gamma1; }border Gamma20(t=0,1) { x=1;
y=t;label=Gamma2; }border Gamma21(t=0,1) { x=1-t; y=1;
label=Gamma2; }int m=10;
mesh Th = buildmesh(Gamma10(m)+Gamma11(m)+Gamma20(m)+Gamma21(m));
plot(Th,wait=true,ps="Th.eps");savemesh(Th,"Th.msh");fespace Vh(Th,P1);
Vh u,v;func f=1;
func g1=0;
func g2=0;
solve Poisson(u,v)=
int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(f*v)-int1d(Th,Gamma2)(g2*v)+on(Gamma1,u=g1); // on(Gamma10,Gamma11,u=g1) とも書ける。
plot(u,wait=1,ps="poisson-kikuchi.eps");//3 次元鳥瞰図//real [int] levels =0.0:0.01:1.0;//plot(u,dim=3,viso=levels,fill=true,wait=true);






curl -O http://nalab.mind.meiji.ac.jp/~mk/program/fem/poisson-kikuchi.edp
FreeFem++ poisson-kikuchi.edp


§2.2 のサンプル・プログラムを叩き台にする。

図 5: poisson-kikuchi.edp — 要素分割と解の等高線
プログラムで生成される、三角形要素分割のデータを出力できるが、それは次のようなものである (しばしば、数値実験の三角形要素分割を、FreeFem++ の出力を拝借している研究者
がいる)。

11


m = 2 のときの “Th.msh”



9 8 8
0 0 1
0 1 2
0 0.5 1
0.5 0 1
1 0 2
0.499999999488 0.499999999488 0
0.5 1 2
1 0.5 21 1 2
7 8 9 01 4 3 04 5 8 06 8 7 04 6 3 04 8 6 02 3 7 07 3 6 0
9 7 2
7 2 2
5 8 2
8 9 2
1 4 1
4 5 1
2 3 1
3 1 1




図 6: m = 2 の時の三角形分割の様子 (Th.msh に対応)
マニュアルで確認していないが、多分次のようなフォーマットになっているのであろう。

12


FreeFem++ のメッシュファイルのフォーマット



総節点数 n 総要素数 m 境界に属する辺の数 N
節点 P1 の x,y 座標とラベル (0 は内点)..
.
節点 Pn の x,y 座標とラベル (0 は内点)要素 e1 の 3 節点の節点番号と 0
要素 e2 の 3 節点の節点番号と 0..
.
要素 em の 3 節点の節点番号と 0境界に属する辺 Q1 とラベル..
.
境界に属する辺 QN とラベル

なお、square()17という関数を使うと、分割まで込めて菊地 [1] の例と同じように出来る。

poissno-kikuchi-square.edp18

// poisson-kikuchi-square.edp// http://nalab.mind.meiji.ac.jp/~mk/program/fem/poisson-kikuchi-square.edp// 菊地文雄, 有限要素法概説, サイエンス社
int m=10;
mesh Th=square(m,m);plot(Th,wait=true,ps="Th.eps");savemesh(Th,"Th.msh");fespace Vh(Th,P1);
Vh u,v;func f=1;
func g1=0;
func g2=0;
solve Poisson(u,v)=
int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(f*v)-int1d(Th,2,3)(g2*v)+on(1,4,u=g1);
plot(u,wait=1,ps="poisson-kikuchi-square.eps");//3 次元鳥瞰図//real [int] levels =0.0:0.01:1.0;//plot(u,dim=3,viso=levels,fill=true,wait=true);B

2 次元 Laplacian の固有値問題を FreeFem++ で解く(ここは書きなおす必要があることが判明したけれど、しばらくは実行する時間が取れない…
書いてあることを信用しないように。)
F. Hecht, O. Pironneau, FreeFem++ a software to solve PDE19 に分かりやすい解説がある(リンク切れ…と思ったら、マニュアルに入った…と安心していたらなくなった…と思ったら復活した)。次に掲げるプログラムはそこに載っているものである。

17square は正方形という意味であるが、square(m,n,[a+(b-a)*x,c+(d-c)*y]) とすると、長方形 [a, b]×[c, d]を分割出来る。
19
https://www.ljll.math.upmc.fr/hecht/ftp/CIMPA/CIMPA-Guadeloupe-FF.pdf

13図 7: m = 10 の時の三角形分割の様子LaplacianEigenvalues.edp

verbosity=10 ;mesh Th=square(20,20,[pi*x,pi*y]);fespace Vh(Th,P2);
Vh u1,u2;
real sigma = 20; // value of the shift// OP = A - sigma B; // the shifted matrix
varf op(u1,u2)= int2d(Th)( dx(u1)*dx(u2) + dy(u1)*dy(u2) - sigma* u1*u2 )+ on(1,2,3,4,u1=0) ; // Boundary condition
varf b([u1],[u2]) = int2d(Th)( u1*u2 ); // no Boundary condition
matrix OP= op(Vh,Vh,solver=Crout,factorize=1);// crout solver because the matrix in not positive
matrix B= b(Vh,Vh,solver=CG,eps=1e-20);// important remark:// the boundary condition is make with exact penalisation:// we put 1e30=tgv on the diagonal term to lock the degre of freedom.
// So take dirichlet boundary condition just on a variationnal form// and not on b variationnanl form.
// because we solve w=OP^-1*B*v
int nev=20;//
real[int] ev(nev); //Vh[int] eV(nev);//

number of computed eigen valeu close to sigma
to store the nev eigenvalue
to store the nev eigenvector

int k=EigenValue(OP,B,sym=true,sigma=sigma,value=ev,vector=eV,
tol=1e-10,maxit=0,ncv=0) ;// return the number of computed eigenvaluefor (int i=0;i<k;i++) {u1=eV[i];real gg = int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1));
real mm= int2d(Th)(u1*u1);
cout<<"----"<< i<<""<<ev[i]<<"err="<<int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1) - (ev[i])*u1*u1)<< " --- "<<endl;
plot(eV[i],cmm="Eigen Vector "+i+" valeur =" + ev[i] ,wait=1,value=1);
}


14



Eigen Vector 0 valeur =5.0002

Eigen Vector 5 valeur =13.0039

IsoValue-0.658137-0.588859-0.519582-0.450304-0.381027-0.311749-0.242472-0.173194-0.103916-0.0346388
0.0346388
0.103916
0.1731940.242472
0.311749
0.381027
0.450304
0.519582
0.588859
0.658137

IsoValue-0.774008-0.692534-0.611059-0.529585-0.44811-0.366635-0.285161-0.203686-0.122212-0.0407373
0.0407373
0.122212
0.203686
0.285161
0.366635
0.44811
0.529585
0.611059
0.692534
0.774008

Eigen Vector 19 valeur =34.0492

IsoValue-0.860704-0.776466-0.692228-0.60799-0.523752-0.439513-0.355275-0.271037-0.186799-0.102561-0.018323
0.065915
0.150153
0.234391
0.318629
0.402867
0.487105
0.571343
0.655582
0.73982

なお、8 行目の +on(1,2,3,4,u=0) を削除すると、Dirichlet 境界条件の代わりに、自然境界条件である Neumann 境界条件となる。

− △ u = λu in Ω,

u = 0 (on ∂Ω)の両辺に試験関数 v をかけて積分し、部分積分すると (Green の積分公式を用いて変形すると)
ZZ
Z∇u · ∇v dx dy = λ uv dx dy.
Ω

すなわち

Ω

ZZ

ZZ
(ux vx + uy vy ) dx dy = λΩuv dx dy.
Ω

これを有限要素法で離散化すると、(行列の) 一般化固有値問題と呼ばれる

Au = λBuという形をした方程式が得られる。ここで B は正値対称行列、A は実対称行列である。
「シフト法で解く」と簡単に書いてあるだけで、詳しいことは書いてない。次のようなことかと推測する。一般化固有値に近いと考えられる実数 σ を与えたとき、Au = λBu の両辺から σBu を引くと、(A − σB)u = (λ − σ)Bu.
ゆえに(A − σB)−1 Bu =

1
u.
λ−σ

1
に属する固有ベクトルであることを示す。
λ−σ−1
ゆえに行列 (A − σB) B について冪乗法を実行すれば、絶対値が最大であるような固有値1が得られる。
µ=
λ−σ
細かいことだが、重要な注意がある。ここでは square() を用いてメッシュ分割をしているが、それを用いずに、自分で境界曲線を定義してメッシュ分割をすると、真の固有値が重複(縮重) する場合であっても、数値計算で得られる近似固有値は重複 (縮重) しない。これは自動メッシュ分割で得られる三角形分割は、正方形の二面体群 (合同変換全体のなす群) D4 の対称性を持たないため、近似固有値が重複しないためである。
これは u が、行列 (A − σB)

−1

の、固有値

(2021/1/9 加筆) 説明する必要が生じて、昨日久しぶりにここを見て、抜本的な書き換えはしないけれど (時間がない)、少し補足することに。
• 行列の固有値問題について、“一般化固有値問題” というキーワードを知っている必要がある。また有限要素法により Au = λBu, A = ((∇ϕi , ∇ϕj )), B = ((ϕi , ϕj )) という行列の一般化固有値問題に帰着される、ということも分かっている必要がある。
15

• 行列の固有値問題についてシフト法 (FreeFem++ の Documentation の表現では、theshifted-inverse mathod) について知っている必要がある。
• FreeFem++ の EigenValue() という関数については、まず、現在の Documentation の§§4.7.46 EigenValue を見ると良い。そこに Parameters (引数と訳すべきか) の説明がある。ただし、maxit, ncv については説明がない。“Other parameters are available formore advanced use: see the ARPACK documentation” とあるので、固有値計算ライブ
ラリィ ARPACK を調べに行く必要がある (ちょっと面倒…昨日はサボった)。

16

17

C

FreeFem++ のマニュアル §3.9 から: Stokes 方程式の cavity ﬂow/****************************************************************************//* This file is part of FreeFEM.
*//**//* FreeFEM is free software: you can redistribute it and/or modify*//* it under the terms of the GNU Lesser General Public License as*//* published by the Free Software Foundation, either version 3 of*//* the License, or (at your option) any later version.
*//**//* FreeFEM is distributed in the hope that it will be useful,*//* but WITHOUT ANY WARRANTY; without even the implied warranty of*//* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*//* GNU Lesser General Public License for more details.
*//**//* You should have received a copy of the GNU Lesser General Public License *//* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>.
*//****************************************************************************/// Parameters
int n = 3; // mesh quality// Meshmesh Th = square(10*n, 10*n);// Fespacefespace Uh(Th, P1b);
Uh u, v;
Uh uu, vv;
fespace Ph(Th, P1);
Ph p, pp;// Problem
solve stokes([u, v, p], [uu, vv, pp])
= int2d(Th)(dx(u)*dx(uu) + dy(u)*dy(uu)+ dx(v)*dx(vv) + dy(v)*dy(vv)+ dx(p)*uu + dy(p)*vv+ pp*(dx(u) + dy(v))-1e-10*p*pp)+ on(1, 2, 4, u=0, v=0)
+ on(3, u=1, v=0);// Plotplot([u,v],p,wait=1);

18

version 4.9 では、/usr/local/ff++/share/FreeFEM/4.9/examples/examples/stokes.edpに置いてある。


cp /usr/local/ff++/share/FreeFEM/4.9/examples/examples/stokes.edp .
FreeFem++ stokes.edp




!
uh =

u1
u2

←

!
u
,
v

!
vh =

v1
v2

←

!
uu
,
vv

ph ← p,qh ← pp,x1 ← x,x2 ← yとすると、
∂u1 ∂v1 ∂u1 ∂v1 ∂u2 ∂v2 ∂u2 ∂v2
+
+
+
dx dy
∂x1 ∂x1 ∂x2 ∂x2 ∂x1 ∂x1 ∂x2 ∂x2
Ωh



ZZ
ZZ 
∂ph
∂u1 ∂u2∂phv1 +v2 dx dy +++qh
dx dy = 0,∂x1∂x2∂x1 ∂x2
Ωh
Ωh

ZZ

すなわち

a(uh , v h ) + (∇ph , v h ) + (qh , ∇ · uh ) = 0 (∀(v h , qh )).
これは(a(uh , v h ) + (∇ph , v h ) = 0 (∀v h ),(qh , ∇ · uh ) = 0 (∀qh )と同値である。
P1b/P1 を採用した結果は、図 8 となる。

図 8: stokes-p1b.edp の結果

P1/P1 にすると上手く計算できないことは常識とされているが、自分でプログラムを書く
場合、なかなかそこまで試す気にはなれない。しかし FreeFem++ でそれをするのは簡単である (図 9)。
19

図 9: もしも P1/P1 とすると…なんだ、これは？
！

図 10: P2/P1 (分割は粗くした)20

P2/P1 は、良く知られているようにきちんと解ける (図 10)。
他にも、悪い冗談のようだが P2/P2 などが気軽に試せる。

D

熱方程式に対する後退 Euler 法

内部で熱が発生する (or 熱が吸収される) 場合の熱方程式の初期値境界値問題(5a)ut (x, t) = △ u(x, t) + f (x) ((x, t) ∈ Ω × (0, ∞)),

(5b)u(x, t) = g1 (x) ((x, t) ∈ Γ1 × (0, ∞)),∂u(x, t) = g2 (x) ((x, t) ∈ Γ2 × (0, ∞)),∂nu(x, 0) = u0 (x) (x ∈ Ω)

(5c)
(5d)

を考える。
変数 t については (後退) 差分近似する。すなわち tn := n∆t, un := u(·, tn ) とおいて、∂uun − un−1(·, tn ) ≒∂t∆tという近似を採用する。
弱形式は


すなわち


un − un−1, v + ⟨un , v⟩ − (f, v) − [g2 , v] = 0.
∆tun+1 , v − (un , v) + ∆t⟨un+1 , v⟩ − ∆t(f, v) − ∆t[g2 , v] = 0.

既に un が “分かっている” とき、これは un+1 についての Poisson 方程式もどきに対する弱形式であり、これまでとほぼ同様にして解くことが出来る。

21


heatB.edp20

// heatB.edp// http://nalab.mind.meiji.ac.jp/~mk/program/fem/heatB.edp// 菊地文雄, 有限要素法概説, サイエンス社の Poisson 方程式の問題の非定常版// 時刻については後退 Euler 法 (陰解法)
int i,m=10;
real Tmax=10, tau=0.01, t;// コメント・アウトすると、m, tau, theta が実行時に変更できるようになる// cout << "m dt theta: ";// cin >> m >> tau >> theta;// cout << "m=" << m << ", tau=" << tau << ", theta=" << theta << endl;
mesh Th=square(m,m);plot(Th,wait=true);fespace Vh(Th,P1);
func f=1;
func g1=0;
func g2=0;
func u0=sin(pi*x)*sin(pi*y);Vh u=u0,uold,v;plot(u,cmm="t=0",wait=1);
problem heat(u,v,init=i)=
int2d(Th)(u*v+tau*(dx(u)*dx(v)+dy(u)*dy(v)))-int2d(Th)(tau*f*v)-int1d(Th,2,3)(tau*g2*v)-int2d(Th)(uold*v)+on(1,4,u=g1);
for (i=0;i<Tmax/tau;i++) {uold=u;t=(i+1)*tau;heat;plot(u,cmm="t="+t,wait=0); // ps="heat"+i+".ps" として保存することも可能}
plot(u,wait=1);






curl -O http://nalab.mind.meiji.ac.jp/~mk/program/fem/heatB.edp
FreeFem++ heatB.edp


ループの制御変数を i という変数にして、problem に ,init=i と書き足すのがミソ。最初は i が 0 であるので、init が False であるが、それ以降は i ̸= 0 であるので、init が True
である、つまり LU 分解済みだ、と指示するわけである。
時間が経過すると、定常解に収束する (Poisson 方程式の解に近く) 様子が見える。

22