2017-11-08

【数値計算】Gauss-Seidel法（ガウス・ザイデル法）のC++コード

数値計算

現在「【数値計算】Gauss-Seidel法（ガウス・ザイデル法）の説明　実践編」を準備中ですが、それに先立ってGauss-Seidel法のC++コードを公開します。私はコーディングの専門家ではないのであまり綺麗ではありませんが悪しからず。皆様の勉強に役立てて頂ければ幸いです。解いている問題は

$\left( \begin{array}{ccc} 3 & 2 & 1 \\ 1 & 4 & 1 \\ 2 & 2 & 5 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 10 \\ 12 \\ 21 \end{array} \right)$
で、答えは
$\left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 1 \\ 2 \\ 3 \end{array} \right)$
となります。実際に動かして遊んでみて下さい。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <vector>

using namespace std;

//inner product//
inline double ip(vector<double> &a, vector<double> &b, int &size)
{
	double value=0.;
	int i;
	for(i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

//Gauss-Seidel
inline void GS(vector<vector<double> > &A,vector<double> &x,vector<double> &b, int &size, double &eps)
{
	int i,k;
	double in,sum,max;
	
	for(k=0;k<1000000;k++)
	{
		max=0.0;
		for(i=0;i<size;i++)
		{
			in=x[i];
			x[i]=(b[i]-(ip(A[i],x,size)-A[i][i]*x[i]))/A[i][i];
			if(max<abs(in-x[i])) max=abs(in-x[i]);
		}
		if(eps>max) break;
	}
}

int main()
{
	int size=3;
	int i;
	double eps=pow(2.0,-50);
	double in,max;
	vector<double> b(size,0.);
	vector<double> x(size,0.);
	vector<vector<double> > A(size,(vector<double>(size,0.)));
	
	A[0][0]=3;A[0][1]=2;A[0][2]=1;
	A[1][0]=1;A[1][1]=4;A[1][2]=1;
	A[2][0]=2;A[2][1]=2;A[2][2]=5;
	b[0]=10;b[1]=12;b[2]=21;
	
	GS(A,x,b,size,eps);

	for(i=0;i<size;i++) cout<<x[i]<<endl;

        return 0;
}

2017-11-01

記事を書いていく順番について

その他

ある事柄について書こうとすると、それを説明するための前提知識についても書かなくてはならなくなって、さらにその前提知識を説明するための前提知識についても書かなくてはならなくなって…、ときりがないのでこれからはどんどん書いていくことにします。もちろん、後で前提知識に関する記事を書いた場合はリンクを張る等していきます。あんまりしっかりやろうとして何も書けない状態よりは、多少がたがたしていてもどんどん書いていきます。

2017-09-24

【数値計算】Gauss-Seidel法（ガウス・ザイデル法）の説明　理論編

数値計算

Gauss-Seidel法とは

連立一次方程式を解くための基本的な反復法です。反復法とは一回の計算で解を求めるのではなく、ある決まった手順を何度も繰り返しながら真の解へと近づいていく方法です。数値計算では有限桁までしか表現できないので、得られる解は誤差を含んだものになります。ちなみに、何故連立一次方程式を解くのかというと、有限要素法などで方程式を離散化すると最終的に連立一次方程式に帰着されるからです。

反復手順

まず解こうとする連立一次方程式を

$Ax=b$

としましょう。ここで、 $A$ は $n \times n$ の行列、 $x$ と $b$ は $n$ 次元のベクトルです。 $A$ と $b$ が既知で、 $x$ が未知です。例えば、 $3 \times 3$ だと

$\left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} b_1 \\ b_2 \\ b_3 \end{array} \right)$

です。この行列 $A$ を

$A=D+L+U$

と分解します。ここで、 $D$ は対角行列（対角成分以外は成分が0の行列）、 $L$ は下三角行列（対角線の下にだけ非零の成分が入っている行列）、 $U$ は上三角行列（対角線の上にだけ非零の成分が入っている行列）です。いかめしく書いていますが、具体例を見ればすぐ分かると思います。例えば、

$A=\begin{pmatrix} 10 & 1 & 1 \\ 2 & 10 & 1 \\2 & 2 & 10 \\ \end{pmatrix}$

としましょう。このとき、 $A$ は、

$A=\begin{pmatrix} 10 & 1 & 1 \\ 2 & 10 & 1 \\ 2 & 2 & 10 \\ \end{pmatrix}=\begin{pmatrix} 10 & 0 & 0 \\ 0 & 10 & 0 \\ 0 & 0 & 10 \\ \end{pmatrix}+\begin{pmatrix} 0 & 0 & 0 \\ 2 & 0 & 0 \\ 2 & 2 & 0 \\ \end{pmatrix}+\begin{pmatrix} 0 & 1 & 1 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ \end{pmatrix}$

のように分解できるよね、と言っているだけです。さて、この分解を最初の連立一次方程式に代入すると

$(D+L+U)x=b$

となります。次に、 $Ux$ を右辺に移項して

$(D+L)x=-Ux+b$

です。最後に、この両辺に $D+L$ の逆行列 $(D+L)^{-1}$ を左からかけて

$x=-(D+L)^{-1}Ux+(D+L)^{-1}b$

を得ます。これをもとにGauss-Seidel法の反復式を作ります。まず、 $k$ ステップ目での $x$ の近似解を $x^k$ としましょう。次のステップ（ $k+1$ ステップ目）での $x$ の近似解を先ほどの式を使って

$x^{k+1}=-(D+L)^{-1}Ux^k+(D+L)^{-1}b$

のように決めることにします。これがGauss-Seidel法の反復式です。適当な初期値 $x^0$ からスタートした数列 $x^0, x^1,\cdots$ は、もし収束するなら $Ax=b$ を満たす $x$ へと収束します。なぜなら、数列 $x^0, x^1,\cdots$ が $x^*$ に収束したとすると

$x^{*}=-(D+L)^{-1}Ux^*+(D+L)^{-1}b$

が成り立ちます。後はこの式を逆に変形していけば

$Ax^*=b$

となるので $x=x^*$ です。

実際の数値計算は有限桁なので、誤差が一定値以下になったら計算を切り上げます。誤差の決め方ですが、自分で決めた誤差の下限 $\epsilon$ （小さい値）に対して

$\max_{1 \leq i \leq n} |x_i^{k+1} - x_i^{k}|\ < \epsilon$

となったら反復を止めることにします。この式の意味は、 $k+1$ ステップ目の近似解と $k$ ステップ目の近似解の差を取ったとき、その差が最大になるものが、あらかじめ定めておいた誤差の下限を下回ったら反復を止めるということです。つまり、近似解の値があまり更新されなくなったら反復を止めよう、と言っているだけです。実は他にも止める基準はあるのですがそれはまた今度。

長くなってきたので一度切りますが、これだけではきっと分からないと思うので、具体的な行列に対する例は次回書きます。

注　Gauss-Seidel法の収束については、結構大変なので別の記事で触れます。

参考

英語で学ぶ数値解析

作者: 陳小君,山本哲朗
出版社/メーカー: コロナ社
発売日: 2002/10/01
メディア: 単行本
この商品を含むブログ (1件) を見る

楽しい反復法

作者: 仁木滉,河野敏行
出版社/メーカー: 共立出版
発売日: 1998/01/01
メディア: 単行本
この商品を含むブログ (1件) を見る

2017-09-23

【語学学習】どうやって新しい外国語を読めるようにするか？

語学学習

私の外国語遍歴？

私は今までに色々な言語に手を出して参りました。大学での第二外学国語のアラビア語に始まり、ラテン語、ドイツ語、デンマーク語、オランダ語、古典ギリシャ語、フランス語、チェコ語、ロシア語、スペイン語、イタリア語。しかし、いろいろやりすぎて何も身に着かない状態が続いていました（それでも様々な言語に触れることが出来てよかったと思っていますが）。そこで、まずはひとつの言語をある程度までしっかり勉強してから他の言語に手を出すことにしました。そこで始めたのがドイツ語です。

文法を身につけましょう

ドイツ語を始めるにあたって、まずは薄い参考書を用いてその言語の全体をおおまかに把握することから始めました。あまり分厚い本に最初から取り組むと挫折しやすいので、最初は薄い参考書から始めましょう。この参考書をとにかく一度読み通します。活用や単語は、こんなのがあるんだ、ぐらいでどんどん進んでいきます。読み終えたらすぐに二週目を始めます。二回目は単語や活用をちゃんと覚えようと努力します。これで基礎文法は大丈夫だと思います。今まで私が語学学習を挫折してきたのはこの二週目の復習が無かったから、だと思います。一回読んだくらいだと人間どんどん忘れてしまいます。欲を言えば三週すれば万全でしょう。

単語もちゃんと覚えましょう

実は、今まで私が語学学習を挫折してきた理由はもう一つあるんです。それは、基礎的な単語を1000語以上覚える努力を完全に怠っていたのです。入門書を一周してだいたいのその言語の特徴を把握するだけで満足していたのです。それではいけません！単語集等を使って覚えましょう。千野栄一先生も著書の中で述べられているように、「語学で必要なのは語彙と文法」なのです。さらに、「基礎語を1000覚えてしまえば初級は卒業でその言語は無に戻ることはない」ということなのです。だから新しい外国語をやる時は、薄い文法書1冊と基礎語1000語を覚えることができる環境を作り出さないといけません。これは例えば、よいテキストや単語集、外国語に取り組む時間、そして気力ないしやる気です。

私のドイツ語は、文法書をやった後、"Mastering German Vocabulary"や『ドイツ基本語5000辞典』を使ってちゃんと単語を覚えたので（後者は10月に終わる予定）、辞書をひきひきWikipediaを読めるようになりました。これで初級は卒業です。初級を卒業したら他の外国語に手を出しても大丈夫です。今、私もフランス語をはじめています。

薄い参考書が終わったら、同じようなレベルの参考書を何冊か読む、のもオススメです。文法の復習にもなるし、場合によってはある本には載っていない文法事項を学習できたりします（入門書は意外とレベルがまちまちでどこまで載っている範囲が全然違う）。さらにここで未知語を拾っておくと1000語の壁を超えるのがかなり楽になります。

外国語上達法 (岩波新書黄版 329)

作者: 千野栄一
出版社/メーカー: 岩波書店
発売日: 1986/01/20
メディア: 新書
購入: 32人クリック: 345回
この商品を含むブログ (83件) を見る

Mastering German Vocabulary: A Thematic Approach (Barron's Vocabulary Series)

作者: Veronika Schnorr
出版社/メーカー: Barrons Educational Series
発売日: 1995/08/01
メディア: ペーパーバック
クリック: 4回
この商品を含むブログ (3件) を見る

ドイツ基本語5000辞典

作者: 岩崎英二郎
出版社/メーカー: 白水社
発売日: 1971/04
メディア: 単行本
この商品を含むブログを見る

そして終わりのない中・上級へ

難しいのが中級からです。これはもうどんどん原文を読んでいくしかありません。なるべく自分の興味がある内容がよいです。Wikipediaはこのためにとても役に立ちます。あらかじめ日本語や英語でその記事を読んでおけば、ドイツ語の記事の内容が想像できます。類推もかなり効くようになります。しかも、専門性の高い記事ほど英語を訳しただけだったりするのでかなり読めます。さらに関連事項も読んでいくとよいです。Wikipediaを読みつつ、未知語を調べ暗記し、わからない文法事項は詳しめの文法書で調べて解決しましょう。やはり、読んでいて自然と出会う語彙がその人にとって本当に必要な語彙である、と思います。あと多少メジャーな言語だと中級向けの解釈の参考書や読本があるのでそれで学習しましょう。基本的に英語以外の言語の参考書はいつ絶版になってAmazonで高額で購入しないといけなくなるかわからないので積極的に入手していきましょう。私も将来やるであろう外国語の参考書を普段から集めています。趣味みたいなもんです。

注　新しい外国語は同時並行で2つ以上勉強しないほうがいいです。経験上どちらもポシャります。

2017-09-22

【有限要素法】有限要素法のプログラミングを勉強するときにどの順序で学ぶのがよいか？【導入編】

数値計算有限要素法

有限要素法のプログラミングを勉強したいと思っている方はたくさんいると思います。私も最初はちんぷんかんぷんでした。Poisson方程式やら移流拡散方程式やらたくさん出てきます。有限要素法や数値計算の本を山ほど仕入れてきてうんうん唸っていました。そんな悪戦苦闘を通して、最近有限要素法がだんだんわかってきました。有限要素法を効率よく勉強するためには、以下の問題に対するプログラムを作成していくのが一番効率がよいと思います。

1次元Poisson方程式（1次元Laplace方程式でも可）
↓
1次元定常移流拡散方程式
↓
1次元非定常拡散方程式
↓
1次元非定常移流拡散方程式
↓
1次元非定常Burgers方程式

まずは1次元から始めるのがよいと思います。もちろん有限要素法が威力を発揮し始めるのは2次元以上ですがシンプルなものから勉強するべきです。

1次元Poisson方程式

そこで最初に1次元Poisson方程式を解くコードを作ります。これは拡散項のある定常問題を解けるようにするためです。ソース項のないLaplace方程式でも可です。ここで大規模な連立一次方程式の解法も学びます。私はGauss-Seidelを愛用していますが、Gaussの消去法も押さえておきましょう。TDMAでもokです。前者が反復法の代表で、後者が直接法の代表です。

1次元定常移流拡散方程式

この問題で拡散項と移流項が存在する定常問題の解き方を学びます。

1次元非定常拡散方程式

この問題では、拡散項だけが存在する非定常問題の解き方を学びます。陽解法、陰解法ともに作ります（ $\theta$ 法で離散化するとよいです）。

1次元非定常移流拡散方程式

次に非定常拡散方程式に移流項を足した問題の解き方を学びます。陽解法、陰解法ともに作ります。これが出来たら線型はもう大丈夫です。非線型問題を解くことは、線型問題を繰り返し解いていくのと同じなので線型問題が解けることは非常に重要です。

1次元非定常Burgers方程式

最後が1次元非定常Burgers方程式です。この問題では、非線型問題の解き方を学びます。Burgers方程式の場合は移流項が非線型項になります。将来Navier-Stokes方程式を解くときにも役立ちます。陽解法、陰解法ともに作ります。陰解法の場合は適当な反復法（Newton法とかPicard反復とか）が必要なのでそれも学びます。

導入なのでざっくりとしか書いていませんが流れはこんな感じです。このあとは風上化編に続きます（執筆中）。コードを書く際に参考になる文献の紹介や各問題の詳細な解説はまた別の記事でやります。

注　1次元移流方程式はリストに入れていませんが、これは風上化を導入した後に解くべきなのであとまわしになっています。

2017-09-20

【語学学習】Wikipediaを使った中・上級者向けの英語学習法

英語語学学習

中・上級者にとって一番大事なことはどんどん原文を読み進めていくことと単語をどんどん覚えることです。ただ単語を覚えるといっても無味乾燥な例文で単語を覚えるのが苦手な方も多いと思います。さらに、中級のうちは英語ならば単語集もありますが、上級になってくるとちょうどよい単語集というのは見当たらなくなってきます。そこで原文（例えばWikipedia英語版の自分の興味のある記事）を読みつつ、出てくる未知語を暗記していくと、非常に効率よく学習を進められます。未知語の暗記にはスマホアプリのAnkiを使います。

http://mathlang.hatenablog.com/entry/2017/09/18/202602

では早速やり方を解説していきます。例えばWikipedia英語版の森鴎外のページを読んでいて知らない単語"acrimonious"が

He divorced her the following year under acrimonious circumstances that irreparably ended his friendship with Nishi.

このように出てきたとします。まずは、Weblioで未知語の意味を調べます。次に未知語と用例、そしてその意味を以下の写真のようにAnkiに打ち込みます。用例はWeblioのものを使ってもよいですし、私は記憶残りやすさから実際の読んでいた文章を用例として使っています。私も毎日この方法でわからない単語、熟語をどんどんAnkiに放り込んでいきます。ノートを使って暗記していたころは、未知語が出てきてもどうせ覚えきれないだろうと考えてやる気がなくなっていましたが、今はAnkiに放り込みさえすればいつか覚えられるので英文をどんどん読み進めることができます！

Wikipediaだけでなく、TimeやBBC Newsを読むのもおすすめです。なかなか高級な語彙が出てきて楽しいです。

ちなみにAnkiで復習する際、私はどうしても単語だけでは思い出せない時だけ用例を読むことにしています。こうすることにより、スピーディな復習ができます。

さらに、この方法は英語などと違って、ちゃんとした単語集のないようなマイナー言語を勉強するのにも適しています。Wikipediaならだいたいのマイナー言語をカバーしています。

2017-08-11

【語学学習】基本語5000辞典

語学学習ドイツ語フランス語ロシア語

ここ二年ほど、毎朝ドイツ語単語帳『ドイツ語基本語5000辞典』を1ページづつ覚えています。この白水社の『基本語5000辞典』にはドイツ語の他に、フランス語、ロシア語、スペイン語とあるのですが、ほぼ全ての単語に例文とその和訳がつけられているという素晴らしい代物です。私もスペイン語以外持っています。もう絶版ですので見つけたら必ず押さえてください。絶対に役に立ちます。私はドイツ語とフランス語は神保町で、ロシア語はAmazonで押さえました。いまだに時々見かけます。飲み会一回パスすれば買えます！

ドイツ基本語5000辞典