2018-09-02

【ドイツ語単語】verschaffen

verschaffen：与える
Thomas verschafft dem Freund eine Wohnung und einen Job.

2018-08-28

【連立一次方程式】双共役勾配法（BiCG法）　C++コード付き

今回は正定値対称行列にしか適用できない共役勾配法（Conjugate Gradient法）をパワーアップした手法である双共役勾配法（Bi-Conjugate Gradient法）のC++コードを公開します。CG法についてはこちらをご覧下さい。

しくみはおいおい説明するとして、とりあえずコードをおいておきます。というより現在進行形で勉強中です！お待ちください。

今回解く連立一次方程式は

$Ax=b$

とあらわされます。ここで、 $A$ は $3 \times 3$ の行列、 $x$ と $b$ は $3$ 次元のベクトルです。 $A$ と $b$ が既知で、 $x$ が未知です。具体的には

$\left( \begin{array}{ccc} 3 & 1 & 1 \\ 1 & 4 & 1 \\ 1 & 1 & 5 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 5 \\ 6 \\ 7 \end{array} \right)$

であり、これを解くと $x_1=1, x_2=1, x_3=1$ となります。コードを書く際に、

数値計算のわざ

作者: 二宮市三
出版社/メーカー: 共立出版
発売日: 2006/02
メディア: 単行本
クリック: 2回
この商品を含むブログ (2件) を見る

の「第4章　大型線型方程式の反復解法」を参考しています。ではコードを以下に示します。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <stdio.h>
#include <time.h>
#include <stdlib.h> 

using namespace std;

int i,j,k,y,x;

//matix vector multiplication Ax=b
inline void mvm(double A[],double x[],double b[], int size)
{
	for(i=0;i<size;i++)
	{
		b[i]=0.;
		for(j=0;j<size;j++)
		{
			b[i]+=A[i*size+j]*x[j];
		}
	}
}

//matix vector multiplication Atx=b (transpose)
inline void tmvm(double At[],double x[],double b[], int size)
{
	for(j=0;j<size;j++)
	{
		b[j]=0.;
		for(i=0;i<size;i++)
		{
			b[j]+=At[j*size+i]*x[i];
		}
	}
}

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

main()
{
	int size=3;
	double alpha,beta;
	double r0,rk,rk1;
	double eps=pow(2.0,-50);
	
	double *r=new double[size];
	double *rs=new double[size];
	double *p=new double[size];
	double *b=new double[size];
	double *ps=new double[size];
	double *Atps=new double[size];
	double *x=new double[size];
	double *Ax=new double[size];
	double *Ap=new double[size];
	double *A=new double[size*size];

	A[0*size+0]=3.;A[0*size+1]=1.;A[0*size+2]=1.;
	A[1*size+0]=1.;A[1*size+1]=4.;A[1*size+2]=1.;
	A[2*size+0]=1.;A[2*size+1]=1.;A[2*size+2]=5.;
	
	b[0]=5.;b[1]=6.;b[2]=7.;

	//IC//
	for(i=0;i<size;i++) x[i]=0.;
	
	mvm(A,x,Ax,size);
	for(i=0;i<size;i++) r[i]=b[i]-Ax[i];
	for(i=0;i<size;i++) p[i]=r[i];
	for(i=0;i<size;i++) ps[i]=r[i];
	for(i=0;i<size;i++) rs[i]=r[i];
	
	r0=sqrt(ip(r,r,size));
	
	//main loop//
	for(k=0;k<10000;k++)
	{
		mvm(A,p,Ap,size);
		alpha=ip(rs,r,size)/ip(ps,Ap,size);
		rk=ip(rs,r,size);
		for(i=0;i<size;i++) x[i]=x[i]+alpha*p[i];
		for(i=0;i<size;i++) r[i]=r[i]-alpha*Ap[i];
		
		rk1=sqrt(ip(r,r,size));
		if(rk1/r0<=eps) break;
		
		tmvm(A,ps,Atps,size);
		for(i=0;i<size;i++) rs[i]=rs[i]-alpha*Atps[i];
		beta=ip(rs,r,size)/rk;
		for(i=0;i<size;i++) p[i]=r[i]+beta*p[i];
		for(i=0;i<size;i++) ps[i]=rs[i]+beta*ps[i];
	}
	
	for(i=0;i<size;i++) cout<<x[i]<<endl;
	
	delete [] r,rs,p,b,ps,Atps,x,Ax,Ap,A;

    return 0;
}

Bi-CG法を安定化した、安定化双共役勾配法（Bi-conjugate gradient stabilized法）はこちらをご覧ください。

2018-08-26

【連立一次方程式】安定化双共役勾配法（Bi-CGSTAB法）　C++コード付き

数値計算連立一次方程式

今回は正定値対称行列にしか適用できない共役勾配法（Conjugate Gradient法）をパワーアップした手法である双共役勾配法（Bi-Conjugate Gradient法）を安定化した、安定化双共役勾配法（Bi-conjugate gradient stabilized法）のC++コードを公開します。ややこしいですね。つまり

CG法→Bi-CG法→Bi-CGSTAB法

のような流れです。CG法についてはこちらをご覧下さい。

しくみはおいおい説明するとして、とりあえずコードをおいておきます。というより現在進行形で勉強中です！お待ちください。

今回解く連立一次方程式は

$Ax=b$

とあらわされます。ここで、 $A$ は $3 \times 3$ の行列、 $x$ と $b$ は $3$ 次元のベクトルです。 $A$ と $b$ が既知で、 $x$ が未知です。具体的には

$\left( \begin{array}{ccc} 1 & 2 & 1 \\ 4 & 1 & 1 \\ 2 & 1 & 1 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 5 \\ 7 \\ 5 \end{array} \right)$

であり、これを解くと $x_1=1, x_2=1, x_3=2$ となります。コードを書く際に、

と
BiCGSTAB法: プログラマーになりたい！

を参考にさせて頂いています。ではコードを以下に示します。

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

using namespace std;

int i,j,k;

//matix vector multiplication Ax=b
inline void mvm(double A[],double x[],double b[], int size)
{
	for(i=0;i<size;i++)
	{
		b[i]=0.;
		for(j=0;j<size;j++)
		{
			b[i]+=A[i*size+j]*x[j];
		}
	}
}

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

main()
{
	int size=3;
	double alpha,beta;
	double r0,rk,rk1,w;
	double eps=pow(2.0,-50);
	
	double *r=new double[size];
	double *p=new double[size];
	double *b=new double[size];
	double *x=new double[size];
	double *Ax=new double[size];
	double *Ap=new double[size];
	double *A=new double[size*size];
	double *rr=new double[size];
	double *s=new double[size];
	double *As=new double[size];
	
	A[0*size+0]=1.;A[0*size+1]=2.;A[0*size+2]=1.;
	A[1*size+0]=4.;A[1*size+1]=1.;A[1*size+2]=1.;
	A[2*size+0]=2.;A[2*size+1]=1.;A[2*size+2]=1.;
	
	b[0]=5.;b[1]=7.;b[2]=5.;
	
	//Initial condition//
	for(i=0;i<size;i++) x[i]=0.;
	
	mvm(A,x,Ax,size);
	for(i=0;i<size;i++) r[i]=b[i]-Ax[i];
	for(i=0;i<size;i++) p[i]=r[i];
	for(i=0;i<size;i++) rr[i]=r[i];
	
	r0=sqrt(ip(r,r,size));
	
	//main loop//
	for(k=0;k<100000;k++)
	{
		mvm(A,p,Ap,size);
		alpha=ip(r,rr,size)/ip(rr,Ap,size);
		rk=ip(r,rr,size);
		
		for(i=0;i<size;i++) s[i]=r[i]-alpha*Ap[i];
		mvm(A,s,As,size);
		w=ip(As,s,size)/ip(As,As,size);
		
		for(i=0;i<size;i++) x[i]=x[i]+alpha*p[i]+w*s[i];
		for(i=0;i<size;i++) r[i]=s[i]-w*As[i];
		
		rk1=sqrt(ip(r,r,size));
		if(rk1/r0<=eps) break;
		
		beta=ip(r,rr,size)/rk*alpha/w;
		for(i=0;i<size;i++) p[i]=r[i]+beta*(p[i]-w*Ap[i]);
	}

	for(i=0;i<size;i++) cout<<x[i]<<endl;
	
	delete [] r,p,b,x,Ax,Ap,A,rr,s,As;

    return 0;
}

2018-08-26

【連立一次方程式】共役勾配法（CG法）　C++コード付き

数値計算連立一次方程式

今回は正定値対称行列に対して適用可能な、共役勾配法（Conjugate Gradient法）のC++コードを公開します。CG法は連立一次方程式の解法のひとつであり、有限回数の反復で解へと収束する面白い方法です（ただし丸め誤差に弱く、必ずしも理論通りにはいってくれない）。連立一次方程式を同値な最適化問題へと置き換えて解を探索していきます。理論は

http://nalab.mind.meiji.ac.jp/~mk/labo/text/linear-eq-2.pdf

http://www.orsj.or.jp/~archive/pdf/bul/Vol.32_06_363.pdf

https://na-inet.jp/nasoft/chap10.pdf

これなら分かる最適化数学―基礎原理から計算手法まで

作者: 金谷健一
出版社/メーカー: 共立出版
発売日: 2005/09/01
メディア: 単行本
購入: 29人クリック: 424回
この商品を含むブログ (42件) を見る

がわかりやすいです。肝は「共役方向」とは何か理解することでしょう。私の理解では楕円に対して円になるような座標変換を施して、円の中心へと向かう方向が「共役方向」です。金谷さんの本と矢部さんのPDFが説明してくれています。

しくみはおいおい説明するとして、とりあえずコードをおいておきます。私も自身のCG法の解説を書く予定です（こうやって書く予定のものがどんどんたまっていきます）。

今回解く連立一次方程式は

$Ax=b$

とあらわされます。ここで、 $A$ は $3 \times 3$ の正定値対称行列、 $x$ と $b$ は $3$ 次元のベクトルです。 $A$ と $b$ が既知で、 $x$ が未知です。具体的には

$\left( \begin{array}{ccc} 1 & 2 & -3 \\ 2 & 5 & -4 \\ -3 & -4 & 8 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} -4 \\ -3 \\ 12 \end{array} \right)$

であり、これを解くと $x_1=0, x_2=1, x_3=2$ となります。ではコードを以下に示します。

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

using namespace std;

int i,j,k;

//matix vector multiplication Ax=b
inline void mvm(double A[],double x[],double b[], int size)
{
	for(i=0;i<size;i++)
	{
		b[i]=0.;
		for(j=0;j<size;j++)
		{
			b[i]+=A[i*size+j]*x[j];
		}
	}
}

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

main()
{
	int size=3;
	double alpha,beta;
	double r0,rk,rk1;
	double eps=pow(2.0,-50);
	
	double *r=new double[size];
	double *p=new double[size];
	double *b=new double[size];
	double *x=new double[size];
	double *Ax=new double[size];
	double *Ap=new double[size];
	double *A=new double[size*size];
	
	
	A[0*size+0]=1.;A[0*size+1]=2.;A[0*size+2]=-3.;
	A[1*size+0]=2.;A[1*size+1]=5.;A[1*size+2]=-4.;
	A[2*size+0]=-3.;A[2*size+1]=-4.;A[2*size+2]=8.;
	
	b[0]=-4;b[1]=-3;b[2]=12;
	
	//Initial condition//
	for(i=0;i<size;i++) x[i]=0.;
	
	mvm(A,x,Ax,size);
	for(i=0;i<size;i++) r[i]=b[i]-Ax[i];
	for(i=0;i<size;i++) p[i]=r[i];
	
	r0=sqrt(ip(r,r,size));
	
	//main loop//
	for(k=0;k<100000;k++)
	{
		mvm(A,p,Ap,size);
		alpha=ip(r,r,size)/ip(p,Ap,size);
		rk=ip(r,r,size);
		
		for(i=0;i<size;i++) x[i]=x[i]+alpha*p[i];
		for(i=0;i<size;i++) r[i]=r[i]-alpha*Ap[i];
		
		rk1=sqrt(ip(r,r,size));
		if(rk1/r0<=eps) break;
		
		beta=ip(r,r,size)/rk;
		for(i=0;i<size;i++) p[i]=r[i]+beta*p[i];
	}

	for(i=0;i<size;i++) cout<<x[i]<<endl;
	
	delete [] r,p,b,x,Ax,Ap,A;

    return 0;
}

CG法を発展させてさらに安定化させた安定化双共役勾配法（Bi-CGSTAB法）についてはこちらをご覧ください。

2018-08-26

【語学学習】語学力（読む力）を中級から上級に伸ばすには？

語学学習英語

ここでは「読み」に関してのみ考えていきます。

中級とは初級文法と初級単語（2,000語以上）を習得した状態としましょう。そして上級をその言語で本や論文を比較的自由に読むことができる（あまり頻繁に辞書や文法書を使わない）状態としましょう。

これはもう多読と精読を合わせた方法をとる一択ですね。つまり丁寧にたくさん文を読む、これだけです。文章を読みながら、文法の確認（中・上級文法を含む）と単語の暗記を繰り返します。そうしているうちに、単語力が増強され、文法は無意識の底へと沈み、半自動化されます。あとはどうやったらこれが可能になるか考えるだけです。私個人の見解としては以下のような方法に至りました。

1. 読む題材としては主にWikipediaから、自分にとって興味があるものを集める。
　日本語から探すのもあり。左のメニューから言語を選択可能。
2. 出てくる未知語は全部覚えること。ただし、日本語でもよくわからない動植物の名前などはその限りではない。
3. 構文が取れない文章が出てきたらGoogle翻訳につっこむ。ただし、英訳にすること。どうやら日本語訳をするよりも精度が高いようである。なので英語ではこの部分は使えない。
4. どうしてもわからない部分は出てくるが、それは頭の隅に残しておいてどんどん文章を読んでいくこと。いつかわかる日が来る。

英語に関してはGoogle翻訳使えませんが、その代わりたくさん参考書が出ているのでそれでカバーすることができると思います。とにかくたくさん読むことです。私もこれでドイツ語とオランダ語を勉強しています。

2018-08-24

【ドイツ語単語】ersetzen

ドイツ語単語ドイツ語語学学習

ersetzen：取り替える、代わりをする
Bewässerung ist die Versorgung des Kulturlandes mit Wasser, um das Wachstum von Pflanzen zu fördern und fehlenden Regen zu ersetzen.

出典
Bewässerung – Wikipedia

2018-08-24

【ドイツ語単語】Bewässerung

ドイツ語単語ドイツ語語学学習

die Bewässerung：灌漑
Bewässerung ist die Versorgung des Kulturlandes mit Wasser, um das Wachstum von Pflanzen zu fördern und fehlenden Regen zu ersetzen.

出典
Bewässerung – Wikipedia