2018-01-16

【英語】"How Starbucks Saved My Life"の一周目【1冊目】

英語語学学習

Michael Gill著、"How Starbucks Saved My Life"の一周目読み終わりました。なかなか未知語が多く、きちんとやると二周目は大変かもしれません。場合によっては二周目は気になる単語だけを拾い、三周目で全て拾うのが正解かもしれません。

How Starbucks Saved My Life: A Son of Privilege Learns to Live Like Everyone Else (English Edition)

作者: Michael Gill
出版社/メーカー: Avery
発売日: 2007/09/20
メディア: Kindle版
この商品を含むブログを見る

翻訳版

ラテに感謝! How Starbucks Saved My Life―転落エリートの私を救った世界最高の仕事

作者: マイケル・ゲイツ・ギル,月沢李歌子
出版社/メーカー: ダイヤモンド社
発売日: 2010/03/05
メディア: 単行本
購入: 4人クリック: 150回
この商品を含むブログ (16件) を見る

2018-01-16

【語学学習】英語とドイツ語とフランス語と

語学学習英語ドイツ語フランス語

英語は辞書なしで小説をスラスラ読む、ドイツ語とフランス語は簡単な文章は辞書無しで、小説は辞書を引き引き読めるようになるのが、今年の目標です。

英語のほうはかなり良いところまで来ているような気がします。ただ、やはり小説を読んでいると見知らぬ単語が現れます。これをできる限り減らしていきたいです。Weblioの語彙力診断をやるとだいたい15,000語から20,000語の間なので、これを30,000語まで増やしたいところです。そのためにとにかく英語の本を読んでいきます。一周目は辞書を引かずに読み、二周目はわからない単語を拾いながら読んでいきます。こうすることによって、多読と精読の利点をどちらも生かすことができます。目標としては100冊ほど読みたいです。あまり厳しくやると続かないので、二周目の読書と対訳本も一冊としてカウントすることとしましょう。

ドイツ語とフランス語は基礎文法は一通りやり、基礎的な単語もまあまあ身に付き始めている、といった状態です。解釈系の教材を用いて、基礎文法、単語の定着を図っていきます。並行して興味のある本も読み進めていければ万々歳です。こちらも最終的に100冊程度読書すれば使える程度に仕上がるのではなかと期待しています。

さらに所々、ラテン語とロシア語も入れていきたいですが、あくまでも優先すべきは上記三言語です。英独仏を固めたいです。

2017-11-30

【数値計算】二次精度風上差分の導出

数値計算差分法

二次精度風上差分の導出をやります。備忘録みたいなものです。点 $i$ を中心に上流側二点の値を使うのが二次精度風上差分です。具体的には、流速が $V_i \geq 0$ のときは $i$ と $i-1$ と $i-2$ の情報を使います。一方、流速が $V_i < 0$ のときは $i$ と $i+1$ と $i+2$ の情報を使います。まず、 $V_i \geq 0$ のときから考えていきます。まず $u_{i-1}$ と $u_{i-2}$ の値を $u_i$ を中心としてTaylor展開します。

$\displaystyle u_{i-1} = u_{i} - \Delta x \frac{u'_{i}}{1!} + (\Delta x)^2 \frac{u''_{i}}{2!} + O[(\Delta x)^3]$
$\displaystyle u_{i-2} = u_{i} - 2\Delta x \frac{u'_{i}}{1!} + (2\Delta x)^2 \frac{u''_{i}}{2!} + O[(\Delta x)^3]$

二階微分の項を消すために $4u_{i-1}$ から $u_{i-2}$ を引きます。すると

$\displaystyle 4u_{i-1} - u_{i-2} = 3u_{i} - 2\Delta x u'_{i} + O[(\Delta x)^3]$

となります。これを $u'_i$ について解くと

$\displaystyle u'_{i} = \frac{3u_{i} - 4u_{i-1} + u_{i-2}}{2\Delta x} + \frac{O[(\Delta x)^3]}{\Delta x} = \frac{3u_{i} - 4u_{i-1} + u_{i-2}}{2\Delta x} + O[(\Delta x)^2]$

となります。つまり二次精度になっています。

同様にして次は、 $V_i < 0$ のときを考えていきます。まず $u_{i+1}$ と $u_{i+2}$ の値を $u_i$ を中心としてTaylor展開します。

$\displaystyle u_{i+1} = u_{i} + \Delta x \frac{u'_{i}}{1!} + (\Delta x)^2 \frac{u''_{i}}{2!} + O[(\Delta x)^3]$
$\displaystyle u_{i+2} = u_{i} + 2\Delta x \frac{u'_{i}}{1!} + (2\Delta x)^2 \frac{u''_{i}}{2!} + O[(\Delta x)^3]$

二階微分の項を消すために $4u_{i+1}$ から $u_{i+2}$ を引きます。すると

$\displaystyle 4u_{i+1} - u_{i+2} = 3u_{i} + 2\Delta x u'_{i} + O[(\Delta x)^3]$

となります。これを $u'_i$ について解くと

$\displaystyle u'_{i} = \frac{-3u_{i} + 4u_{i+1} - u_{i+2}}{2\Delta x} + \frac{O[(\Delta x)^3]}{\Delta x} = \frac{-3u_{i} + 4u_{i+1} - u_{i+2}}{2\Delta x} + O[(\Delta x)^2]$

となります。つまり二次精度になっています。

まとめると

$\displaystyle \left.V\frac{\partial u}{\partial x}\right|_i = \begin{cases} \displaystyle V_i\frac{3u_{i} - 4u_{i-1} + u_{i-2}}{2\Delta x} & (V_i \geq 0)\\ \displaystyle V_i\frac{-3u_{i} + 4u_{i+1} - u_{i+2}}{2\Delta x} & (V_i < 0) \end{cases}$
となります。これが二次精度の風上差分です。

2017-11-28

【数値計算】風上差分の中心差分+数値拡散表示

数値計算差分法

風上差分の中心差分+数値拡散表示の導出をやります。ポイントは風向きの正負によらない表示です。まず、 $u_i$ を節点 $i$ における濃度ないしは熱、 $V_i$ を節点 $i$ における流速、 $x$ を空間座標としたとき、風上差分は

$\displaystyle \left.V\frac{\partial u}{\partial x}\right|_i = \begin{cases} \displaystyle V_i\frac{u_i - u_{i-1}}{\Delta x} & (V_i \geq 0)\\ \displaystyle V_i\frac{u_{i+1} - u_i}{\Delta x} & (V_i < 0) \end{cases}$
のように書けます。風上側（情報が伝わってくる側）の値を使うことで安定的に計算を行うことができます。ここで

$\displaystyle \frac{V_i+|V_i|}{2}= \begin{cases} \displaystyle V_i & (V_i \geq 0)\\ \displaystyle 0 & (V_i < 0) \end{cases}$ 　 $\displaystyle \frac{V_i-|V_i|}{2}= \begin{cases} \displaystyle 0 & (V_i \geq 0)\\ \displaystyle V_i & (V_i < 0) \end{cases}$

というちょっとテクニカルな関係を使います。落ち着いて絶対値をはずしてみればたいしたことはありません。上の関係をよく見ると、 $V_i \geq 0$ のとき $V_i=\frac{V_i+|V_i|}{2}$ に、 $V_i < 0$ のとき $V_i=\frac{V_i-|V_i|}{2}$ になっていることに気づきます（よく見て下さい）。なので、これらを風上差分の式の $V_i$ に代入してみましょう。すると

$\displaystyle \left.V\frac{\partial u}{\partial x}\right|_i = \begin{cases} \displaystyle \frac{V_i+|V_i|}{2} \frac{u_i - u_{i-1}}{\Delta x} & (V_i \geq 0)\\ \displaystyle \frac{V_i-|V_i|}{2} \frac{u_{i+1} - u_i}{\Delta x} & (V_i < 0) \end{cases}$
となります。さらに、 $\boldsymbol{\frac{V_i+|V_i|}{2} \frac{u_i - u_{i-1}}{\Delta x}}$ と $\boldsymbol{\frac{V_i-|V_i|}{2} \frac{u_{i+1} - u_i}{\Delta x}}$ を足してしまいます！するとそれは風上差分と等しくなっているのです！すなわち

$\displaystyle \left.V\frac{\partial u}{\partial x}\right|_i = \frac{V_i+|V_i|}{2} \frac{u_i - u_{i-1}}{\Delta x} + \frac{V_i-|V_i|}{2} \frac{u_{i+1} - u_i}{\Delta x}$
です。何故そんなことができるか？それは場合分けして考えるとわかります。まず $V_i \geq 0$ のときは

$\displaystyle \frac{V_i+|V_i|}{2}=V_i$ 　かつ　 $\displaystyle \frac{V_i-|V_i|}{2}=0$

なので和の二項目が消えて風上差分と等しくなります。また $V_i < 0$ のときは

$\displaystyle \frac{V_i+|V_i|}{2}=0$ 　かつ　 $\displaystyle \frac{V_i-|V_i|}{2}=V_i$

なので和の一項目が消えて風上差分と等しくなります。

さらに式を $V_i$ と $|V_i|$ で整理すると結局

$\displaystyle \left.V\frac{\partial u}{\partial x}\right|_i = V_i \frac{u_{i+1} - u_{i-1}}{2 \Delta x} - \frac{|u_i| \Delta x}{2} \frac{u_{i+1} - 2 u_i + u_{i-1}}{(\Delta x)^2}$
となります。この式の一項目は、一階微分に対する中心差分に、二項目は二階微分に対する中心差分となっているので、風上差分の中心差分+数値拡散表示を得ることができました。これは風向きの正負によらない表示でもあります。

2017-11-26

【語学学習】ラテン語の初級を学ぶためにどんな本を読んだらいいか？

ラテン語語学学習

河島思朗著『基本から学ぶラテン語』を読み終わりました。情報が詰め込まれすぎていない、活用表のレイアウトが非常に見やすい、例文が長すぎず複雑すぎない、初級の文法はきちんと押さえてある、という素晴らしい本でした。ラテン語をはじめようと思っている方に是非おすすめした一冊です。難易度的にはこの本が一番やさしいのではないでしょうか（文法を網羅していない紹介的な本は除いて）。大西英文著『はじめてのラテン語』や山下太郎著『しっかり学ぶ初級ラテン語』よりもとっつきやすく、通読できる可能性が高いのではないでしょうか。有田潤著『初級ラテン語入門』もラテン語の文をどんどん読ませる素晴らしい本なのですが、初級文法を網羅していないのでこれで初級文法を終わりにすることができないのです。

学習順序としては、読み物的な本、初級文法を概観する本、初級文法を網羅した本という順序で

1, 逸身喜一郎著『ラテン語のはなし―通読できるラテン語文法』、
　小林標著『ラテン語の世界―ローマが残した無限の遺産』、
　小倉博行著『ラテン語のしくみ』
2, 河島思朗著『基本から学ぶラテン語』
3, 大西英文著『はじめてのラテン語』
4, 有田潤著『初級ラテン語入門』
5, 山下太郎著『しっかり学ぶ初級ラテン語』

のような感じでしょうか。私は次は『しっかり学ぶ初級ラテン語』をやります。実は一度やって半分で挫折してしまいました。後半の活用、曲用ラッシュに消化不良に陥りやる気がなくなってしまったのです…でも、今回は初級文法は概観してきたので読み切れるはず。頑張ります！『しっかり学ぶ初級ラテン語』を読み切り、消化したら初級文法は一応終わりということにしてよいでしょう。

基本から学ぶラテン語

作者: 河島思朗
出版社/メーカー: ナツメ社
発売日: 2016/07/14
メディア: 単行本（ソフトカバー）
この商品を含むブログ (3件) を見る

はじめてのラテン語 (講談社現代新書)

作者: 大西英文
出版社/メーカー: 講談社
発売日: 1997/04/18
メディア: 新書
購入: 5人クリック: 48回
この商品を含むブログ (36件) を見る

初級ラテン語入門

作者: 有田潤
出版社/メーカー: 白水社
発売日: 1964/01/01
メディア: 単行本
クリック: 2回
この商品を含むブログ (15件) を見る

しっかり学ぶ初級ラテン語 (Basic Language Learning Series)

作者: 山下太郎
出版社/メーカー: ベレ出版
発売日: 2013/08/20
メディア: 単行本
この商品を含むブログ (4件) を見る

ラテン語のはなし―通読できるラテン語文法

作者: 逸身喜一郎
出版社/メーカー: 大修館書店
発売日: 2000/12/01
メディア: 単行本
購入: 4人クリック: 18回
この商品を含むブログ (29件) を見る

ラテン語の世界―ローマが残した無限の遺産 (中公新書)

作者: 小林標
出版社/メーカー: 中央公論新社
発売日: 2006/02/01
メディア: 新書
購入: 4人クリック: 41回
この商品を含むブログ (47件) を見る

ラテン語のしくみ (言葉のしくみ)

作者: 小倉博行
出版社/メーカー: 白水社
発売日: 2014/11/19
メディア: 単行本（ソフトカバー）
この商品を含むブログ (3件) を見る

全てのラテン語の教科書について言えると思うのですが、接続法が出てきてからがしんどいです。何も知らずに読んでいくと、接続法が出てきた途端に覚える活用の量が二倍になって絶望してしまうのです。もっとちょびちょび接続法を出していく方法で書いたらすごい良いラテン語の入門書が書けるかもしれません。

2017-11-24

【数値計算】流束制限関数とは何だろう？そのアイディア

数値計算差分法

流束制限関数というのは、数値流体力学において双曲型の方程式を効率良く解くために考案された手法です。以下にそのアイディアを示します。

移流方程式に対して高次精度の差分を使うと衝撃波（解の勾配がきつくなる所）
近辺で数値解が振動してしまう、つまり正しい解を得られなくなってしまう。
↓
一次精度の風上差分を使うと衝撃波近辺でも数値解は振動しないが、
数値拡散が強く衝撃波がなまってしまう。
↓
高次精度の差分と一次精度の風上差分のよいとこどりをした手法を作ろう。

こうしたアイディアのもと考案されたのが流束制限関数（flux limiter）です。例えばminmodとかsuperbeeとか色々な種類のlimiterがあります。流束制限関数を用いて風上差分と高次の差分をつなぎます。衝撃波の近辺では風上差分で、十分解が滑らかな領域では高次の差分に切り替えて計算を行います。

次回は、数値流束の話、minmodとかsuperbeeの具体的な形やその実装を見ていきましょう。

2017-11-21

【有限要素法】1次元非定常拡散方程式をGalerkin法で解く　C++コード付き

数値計算有限要素法

今回は1次元非定常拡散方程式を有限要素法（Galerkin法）で解いていきます。C++コード付きです。1次元非定常拡散方程式は、前に書いた記事「有限要素法のプログラミングを勉強するときにどの順序で学ぶのがよいか？」で紹介した、有限要素法を学ぶ場合に三番目に解くべき方程式です。

1次元非定常拡散方程式とは

$\displaystyle \frac{\partial u}{\partial t}=D\frac{\partial^2 u}{\partial x^2}, \quad x \in [0, L]$

のような偏微分方程式のことをいいます。ここで、 $u(t,x)$ は熱や溶質の濃度、 $D$ は拡散の効果を表す拡散係数と解釈することができます。熱またはある溶質が時間が経過するにつれてどのような分布をとるか（拡散していくか）求めるのがこの問題です。 $u$ は時間方向と空間方向ともに変化するので二変数関数 $u(t,x)$ となっています。前に説明した1次元Poisson方程式を時間発展させたものにあたります。境界条件は前回と同じく3パターン用意しました。すなわち、

Case 1, $x=0$ と $x=L$ でDirichlet条件を課す
（ $u(t,0)$ と $u(t,L)$ の値を指定）
Case 2, $x=0$ でNeumann条件、 $x=L$ でDirichlet条件を課す
（ $\frac{ \partial u }{ \partial x }(t,0)$ と $u(t,L)$ の値を指定）
Case 3, $x=0$ でDirichlet条件、 $x=L$ でNeumann条件を課す
（ $u(t,0)$ と $\frac{ \partial u }{ \partial x }(t,L)$ の値を指定）

です。それぞれコード上ではbc=1、bc=2、bc=3に対応しています。連立一次方程式の解法はGauss-Seidel法またはTDMAを用いています。好きなほうを選べます。詳しくは以下の記事を参考にして下さい。

空間方向の離散化はGalerkin法、時間方向は $\theta$ 法で離散化しており、今回は $\theta=0.5$ のCrank-Nicolson法を用いています。ちなみに、 $\theta=0$ のとき前進Euler法（陽解法、explicit scheme）に、 $\theta=1$ のときに後退Euler法（完全陰解法、implicit scheme）になります。これらについても今後まとめる予定です。

計算条件については以下のように設定しています。まず、拡散方程式を解く領域は $x\in[0,1$ ]、要素数は100、節点数は101で等分割（ $\Delta x=1/100$ ）、時間刻みは $\Delta t=0.001$ 、計算終了時刻は $t=10$ 、拡散係数は $D=0.01$ 、 $\theta=0.5$ 、境界条件はbc=1（両側Dirichlet条件）で $u(t,0)=0$ と $u(t,1)=0$ 、初期条件は頂点の値が0.5の三角形（下の図参照）です。式で書くと $f(x)=\min(x,1-x)$ です。初期時刻に三角形のような熱の分布を与え、棒の両側をずっと0℃で冷やしたときに、棒の熱はある時間位どのような分布をとるか、という問題です。答えは"Answer_****.txt"（****にはステップ数が入る）に出力されます。1000ステップごとに値が出力されます。一行目が $x$ 座標、二行目がその時刻における $u$ の分布となっています。

コードです。

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

using namespace std;

inline void GS(double AD[],double AL[],double AR[],double x[],double xx[],double b[], int &Node, double &eps)
{
	int i,k;
	double in,max;
	
	for(k=0;k<100000;k++)
	{
		max=0.0;
		for(i=0;i<Node;i++)
		{
			in=xx[i];
			xx[i]=(b[i]-(AL[i]*xx[i-1]+AR[i]*xx[i+1]))/AD[i];
			if(max<abs(in-xx[i])) max=abs(in-xx[i]);
		}
		
		if(eps>max) break;
	}
}

//TDMA//a:diagonal,c:left,b:right,
inline void tdma(double a[], double c[], double b[], double d[], double x[], int size)
{
	int i;
	double *P=new double[size];
	double *Q=new double[size];
	
	//first step//
	P[0]=-b[0]/a[0];
	Q[0]=d[0]/a[0];
	
	//second step//
	for(i=1;i<size;i++)
	{
		P[i]=-b[i]/(a[i]+c[i]*P[i-1]);
		Q[i]=(d[i]-c[i]*Q[i-1])/(a[i]+c[i]*P[i-1]);
	}
	
	//third step, backward//
	x[size-1]=Q[size-1];

	for(i=size-2;i>-1;i=i-1)
	{
		x[i]=P[i]*x[i+1]+Q[i];
	}
	
	delete [] P,Q;
}

inline void mat(double AD[],double AL[],double AR[],double BD[],double BL[],double BR[],double b[],double f[], int &Ele, double &dx, double &Diff, double &dt, double &theta, int &Node)
{
	int i;
	
	//initialization//
	for(i=0;i<Node;i++)
	{
		AD[i]=0.0;AL[i]=0.0;AR[i]=0.0;BD[i]=0.0;BL[i]=0.0;BR[i]=0.0;
	}
	
	for(i=0;i<Ele;i++)
	{
		//A//
		//teporal//
		AD[i]+=dx*2.0/6.0/dt;
		AR[i]+=dx*1.0/6.0/dt;
		AL[i+1]+=dx*1.0/6.0/dt;
		AD[i+1]+=dx*2.0/6.0/dt;
		
		//diffusion//
		AD[i]+=theta*Diff/dx;
		AR[i]+=theta*-Diff/dx;
		AL[i+1]+=theta*-Diff/dx;
		AD[i+1]+=theta*Diff/dx;
		
		//B//
		//temporal//
		BD[i]+=dx*2.0/6.0/dt;
		BR[i]+=dx*1.0/6.0/dt;
		BL[i+1]+=dx*1.0/6.0/dt;
		BD[i+1]+=dx*2.0/6.0/dt;
		
		//diffusion//
		BD[i]+=-(1.0-theta)*Diff/dx;
		BR[i]+=-(1.0-theta)*-Diff/dx;
		BL[i+1]+=-(1.0-theta)*-Diff/dx;
		BD[i+1]+=-(1.0-theta)*Diff/dx;
	}
}

inline void boundary(int &bc,double AD[],double AL[],double AR[],double BD[],double BL[],double BR[],double b[],int &Node)
{
	if(bc==1)
	{
		AD[0]=1.0;AR[0]=0.0;BD[0]=1.0;BR[0]=0.0;
		AL[Node-1]=0.0;AD[Node-1]=1.0;BL[Node-1]=0.0;BD[Node-1]=1.0;
	}
	if(bc==2)
	{
		AL[Node-1]=0.0;AD[Node-1]=1.0;BL[Node-1]=0.0;BD[Node-1]=1.0;
	}
	if(bc==3)
	{
		AD[0]=1.0;AR[0]=0.0;BD[0]=1.0;BR[0]=0.0;
	}
}

inline void text(int &i,double x[], double &dx, int &Node)
{
	int j;
	
	stringstream ss;
	string name;
	ofstream fo;
	
	ss<<i;
	name=ss.str();
	name="Answer_" + name + ".txt";
	fo.open(name.c_str ());
	
	for(j=0;j<Node;j++)
	{
		fo<<dx*float(j)<<" "<<x[j]<<endl;
	}
}

int main()
{
	int i,j;
	int Ele=100;
	int Node=Ele+1;
	double LL=1.0;
	double dx=LL/Ele;
	double dt=0.001;
	double NT=10000;
	double eps=pow(2.0,-50);
	double Diff=0.01;
	double theta=0.5;
	int bc=1;//bc=1 both Diriclet bc=2 left Neumann bc=3 right Neumann
	double D0=0.0;
	double D1=0.0;
	double N0=0.0;
	double N1=0.0;
	
	double *b=new double[Node];
	double *x=new double[Node];
	double *xx=new double[Node];
	double *f=new double[Node];
	double *AD=new double[Node];
	double *AL=new double[Node];
	double *AR=new double[Node];
	double *BD=new double[Node];
	double *BL=new double[Node];
	double *BR=new double[Node];
	
	ofstream fk;
	fk.open("Answer_0.txt");
	
	//initial condition//
	for(i=0;i<Node;i++)
	{
		x[i]=min(dx*float(i),1.0-dx*float(i));
		fk<<dx*float(i)<<" "<<x[i]<<endl;
	}
	
	mat(AD,AL,AR,BD,BL,BR,b,f,Ele,dx,Diff,dt,theta,Node);
	boundary(bc,AD,AL,AR,BD,BL,BR,b,Node);
	
	for(i=1;i<=NT;i++)
	{
		//for b//
		if(bc==1)
		{
			b[0]=D0;
			b[Node-1]=D1;
		}
		if(bc==2)
		{
			b[Node-1]=D1;
			b[0]=BD[0]*x[0]+BR[0]*x[1]-N0*Diff;
		}
		if(bc==3)
		{
			b[0]=D0;
			b[Node-1]=BL[Node-1]*x[Node-2]+BD[Node-1]*x[Node-1]+N1*Diff;
		}
		
		for(j=1;j<Node-1;j++)
		{
			b[j]=BL[j]*x[j-1]+BD[j]*x[j]+BR[j]*x[j+1];
		}

		//GS or TDMA//
		GS(AD,AL,AR,x,xx,b,Node,eps);
		//tdma(AD,AL,AR,b,xx,Node);
		for(j=0;j<Node;j++) x[j]=xx[j];
		
		if(i%1000==0)
		{
			cout<<i<<endl;
			text(i,x,dx,Node);
		}
	}
	
	delete[] b,x,xx,f,AD,AL,AR,BD,BL,BR;
	
	return 0;
}

以下が計算結果の図です。初期に与えた三角形の熱の分布が時間が経過するにつれてどんどんなまっていく様子が観察できます。両端を冷やしているので棒がどんどん冷めていくのです。もっと時間進行させると領域全体で熱の分布が0に収束していくのが観察できます。

f:id:mutsumunemitsutan:20181022213936p:plain:w500

ゆくゆくはこのブログだけで有限要素法やいろいろな語学を学べる、self-contained（必要なものが全てそろった）なものにしていきたいと考えています。なので、記事を読んでいて説明が不足しているなと思われたり、わかりづらいなと思われる箇所がありましたら報せて頂けると幸いです。

離散化の詳細はお待ちください。