C++の入出力ファイル サンプル

#include<iostream>
#include<fstream>
#include<string>
using namespace std;

int main(){
ifstream ifs;
ofstream k;
ofstream q;
int x;
for (int i = 1; i < 4; i++){
ifs.open("test.txt");
k.open("aa" + to_string(i) + ".txt");
for (int j = 0; j < 3; j++){
ifs >> x;
k << i*x << endl;
}
ifs.close();
k.close();
ifs.open("test.txt");
q.open("q" + to_string(i) + ".txt");
for (int j = 0; j < 3; j++){
ifs >> x;
q << i * 2 * x << endl;
}
ifs.close();
q.close();
}
return 0;
}

g++ example.cpp -std=c++11


Makefile example

test:testfile
./testfile
testfile: test.cpp
g++ test.cpp -std=c++11 -o testfile

三次元空間における平均曲率は単位法線ベクトルの発散となることの概要

$$ \bm{p}=(x(u,v),y(u,v),z(u,v)) $$
を考える
$ (u,v)=(a,b) $における接平面は
$$ \bm{p}_u(a,b),\bm{p}_v(a,b) $$が張る平面となるので$(a,b)$における単位法線ベクトルは
$$ \bm{e}=\cfrac{\bm{p}_u(a,b)\times\bm{p}_v(a,b)}{|\bm{p}_u(a,b)\times\bm{p}_v(a,b)|} $$

曲率の定義
弧長変数$s$をもちいて$\bm{p}''(s)$の長さを曲率とする
曲面の法ベクトル方向の$\bm{p}''(s)$の長さを法曲率とする

$ \bm{p}_0 $から出る接平面上のベクトルをとると、ベクトルのとり方は$ \bm{p}_0 $を360度囲むように単位円上にとり方があるので、ベクトルのとり方に寄る法曲率の最大値と最小値の平均を平均曲率とする。接平面上のベクトルを選ぶ上では長さが同じという条件が必要であり、ラグランジュの未定乗数法を用いることが出来る。

$$ \frac{\kappa_1+\kappa_2}{2}=\cfrac{EN+GL-2FM}{2(EG-F^2)}$$
となる。これは
$$ \frac{\kappa_1+\kappa_2}{2}=\cfrac{(p_u\cdot p_u)(p_{vv}\cdot e) +(p_v\cdot p_v)(p_{uu}\cdot e)-2(p_u\cdot p_v)(p_{uv}\cdot e)}{2((p_u\cdot p_u)(p_v\cdot p_v)-(p_u\cdot p_v)^2)}$$

また、$$|\bm{p}_u(a,b)\times\bm{p}_v(a,b)|^2=|\bm{p}_u(a,b)|^2|\bm{p}_v(a,b)|^2\sin^2\theta$$
$$ =|\bm{p}_u(a,b)|^2|\bm{p}_v(a,b)|^2-(\bm{p}_u(a,b)\cdot\bm{p}_v(a,b))^2 $$
を$ X^2 $とすると、
$$ \nabla\cdot \bm{e}=$$

ラグランジュ補間による差分法の2次精度微分

ラグランジュ補間による差分法の2次精度微分


関数$ f $の数値微分$f'(x_1)$を行う方法のメモ．

$$ (x_0,f_0),(x_1,f_1),(x_2,f_2)$$
が分かるとする($x_0<x_1<x_2$, $ f_0,f_1,f_2$は定数)

上記3点を通る2次関数は
$$g= f_0\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}+f_1\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}+f_2\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}$$
である．代入すると自明である．
$$\Delta x_0=x_1-x_0, \Delta x_1=x_2-x_1$$
とおく．
$g$を微分して$x=x_1$を代入すると
$$ f'(x_1)\approx g'(x_1)$$であり，
$$g'(x_1)=-f_0\frac{\Delta x_1}{\Delta x_0(\Delta x_0+\Delta x_1)}+f_1\frac{(\Delta x_1-\Delta x_0)}{\Delta x_0\Delta x_1} +f_2\frac{\Delta x_0}{\Delta x_1(\Delta x_0+\Delta x_1)}$$
\begin{equation}=\cfrac{f_2(\Delta x_0)^2-f_0(\Delta x_1)^2+f_1((\Delta x_1)^2-(\Delta x_0)^2)}{\Delta x_0 \Delta x_1(\Delta x_0+\Delta x_1)}\end{equation}
となる　これが2次精度の公式である．

次に2次精度であることを確認する．
テイラー展開をして
$$f_2=f_1+\Delta x_1(f'(x_1))+\frac{(\Delta x_1)^2}{2!}f''(x_1)+\mathcal{O}((\Delta x_1)^3)$$
$$f_0=f_1-\Delta x_0(f'(x_1))+\frac{(\Delta x_0)^2}{2!}f''(x_1)+\mathcal{O}((\Delta x_0)^3)$$

これらを(1)に代入する
(1)の分子は
$$ f_1(\Delta x_0)^2+\Delta x_1(\Delta x_0)^2(f'(x_1))+\frac{(\Delta x_1)^2 (\Delta x_0)^2}{2!}f''(x_1)+\mathcal{O}((\Delta x_1)^3(\Delta x_0)^2)-f_1(\Delta x_1)^2+\Delta x_0(\Delta x_1)^2(f'(x_1))-\frac{(\Delta x_0)^2(\Delta x_1)^2}{2!}f''(x_1)-\mathcal{O}((\Delta x_0)^3(\Delta x_1)^2)+f_1((\Delta x_1)^2-(\Delta x_0)^2)
$$
$$=\Delta x_0\Delta x_1(\Delta x_0+\Delta x_1)(f'(x_1))+\mathcal{O}\left((\Delta x_0)^2(\Delta x_1)^2(\Delta x_1+\Delta x_0)\right)$$
よって
$$(1)=f'(x_1)+\mathcal{O}\left(\Delta x_0\Delta x_1\right)$$

なので2次精度を満たす！

C++で最小二乗法

|Ax-b|^2を最小にするxは、偏微分で正規方程式をもとめると、
A^tAx=A^tb
の解となる
AをQR分解すると
A^tAx=R^tQ^tQRx=R^tRx
A^tb=R^tQ^tb
となる
(R^t)^-1を両辺にかけて
Rx=Q^tbを求めれば良い
Rは三角なので簡単に求められる

(0, -2)

(1, 2,)

(2, 5)

(3, 11)

(4, 1)に対して3次式で近似してみる

次はQR分解で固有値求める予定
最小二乗法のプログラムは↓

#include<iostream>

using namespace std;

/*　m行n列　*/

double** init_matrix(int m, int n) {

double** A = new double*[m];

for (int i = 0; i < m; i++) {

A[i] = new double[n]();

for (int j = 0; j < n; j++) {

cin >> A[i][j];

}

}

return A;

}

void print_matrix(double** A, int n, int m) {

cout << endl;

cout << endl;

for (int i = 0; i < m; i++) {

for (int j = 0; j < n; j++) {

cout << A[i][j] << " ";

}

cout << endl;

}

}

/*Aの左下三角のそれぞれの列ににハウスホルダー変換のuが入る　y[i]にはA[i][i]がはいり，bで後退代入するときに使う*/

int qr(double** A, int n, int m, double* y) {

double x, z, w;

for (int i = 0; i < n; i++) {

x = 0;

for (int j = i; j < m; j++) {

x += A[j][i] * A[j][i];

}

if (x < 0.00001) {

return 1;

}

x = sqrt(x);

if (A[i][i] > 0) {

y[i] = -x;

}

else {

y[i] = x;

}

z = 1 / sqrt(y[i] * (y[i] - A[i][i]));

/*　ハウスホルダー変換のuをAの左下三角の列に格納　*/

A[i][i] = z*(A[i][i] - y[i]);

for (int j = i + 1; j < m; j++) {

A[j][i] = A[j][i] * z;

}

/* Aの右下にハウスホルダー変換を行う計算

　 H=E-uu^tに対してHv=v-(u,v)uとなる．

   まず，w=(u,v)計算する．

*/

for (int j = i + 1; j < n; j++) {

w = 0;

for (int k = i; k < m; k++) {

w += A[k][j] * A[k][i];

}

for (int k = i; k < m; k++) {

A[k][j] -= w*A[k][i];

}

}

}

return 0;

}

void lsm(double** A,double* y, double* b,int m,int n) {

/* Aの右下にハウスホルダー変換を行う計算

　 H=E-uu^tに対してHb=b-(u,b)uとなる．

  まず，w=(u,v)計算する．

  */

double w;

for (int i = 0; i < n; i++) {

w = 0;

for (int j = i; j < m; j++) {

w += b[j] * A[j][i];

}

for (int j = i; j < m; j++) {

b[j] -= w*A[j][i];

}

}

cout << " b  print" << endl;

for (int i = 0; i < n; i++) {

cout << b[i] << " ";

}

cout << endl;

/*下から代入して解く*/

for (int i = n - 1; i >= 0; i--) {

for (int j = i + 1; j < n; j++) {

b[i] -= b[j] * A[i][j];

}

b[i] = b[i] / y[i];

}

}

int main() {

int n, m;

cin >> m >> n;

double* y=new double[n]();

double** A = init_matrix(m, n);

double* b = new double[m];

for (int i = 0; i < m; i++) {

cin >> b[i];

}

print_matrix(A, n, m);

qr(A, n, m, y);

print_matrix(A, n, m);

cout << "y" << endl;

 for (int i = 0; i < m; i++) {

cout << y[i] << " ";

}

cout << endl;

lsm(A, y, b, m, n);

cout << endl;

cout << "b" << endl;

for (int i = 0; i < n; i++) {

cout<< b[i]<<endl;

}

return 0;

}

結果

1 0 0 0 

1 1 1 1 

1 2 4 8 

1 3 9 27 

1 4 16 64 

matend

1.203 -4.47214 -13.4164 -44.7214 

0.371748 -1.05867 12.6491 48.6991 

0.371748 0.184608 -1.33233 22.4499 

0.371748 0.48331 -0.472748 -1.13192 

0.371748 0.782012 0.0374653 0.847794 

matend

y

-2.23607 3.16228 3.74166 3.79473 2.10535e-314 

 b  print

-7.60263 4.74342 -6.68153 -4.74342 

b

-1.67143

-2.10714

5.71429

-1.25

レイノルズ平均方程式の導出

乱流を扱うときに使うモデルです．
速度を平均と変動成分に分けます．
レイノルズ分解
$$ u_i=\bar{u}_i+u'_i $$
流れが定常→時間平均
一様な流れ→空間平均
一般→アンサンブル平均

アンサンブル平均
$$\bar{u}=\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{i=1}^Nu_i$$
次を満たすアンサブル平均をレイノルズ平均という
$$　\overline{u'}_i=0$$
$$\bar{\bar{u}}_i=\bar{u}_i$$
$$\overline{u'_i\bar{u}_j}=0$$
次が成りたつ
$$ \overline{\frac{\partial u_i}{\partial x_j}}=\frac{\partial \bar{u}_i}{\partial x_j}$$
$$ \overline{\frac{\partial u'_i}{\partial x_j}}=\frac{\partial \bar{u'}_i}{\partial x_j}=0$$
非圧縮流体の連続の式は$$ \frac{\partial u_i}{\partial x_i}=0$$
なのでレイノルズ平均を取ると$$\frac{\partial \overline{u}_i}{\partial x_i}=0$$
両辺引くと$$\frac{\partial u'_i}{\partial x_i}=0$$
ナビエ・ストークス方程式は
$$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho}\frac{\partial p}{\partial x_i}+\nu \frac{\partial^2 u_i}{\partial x_j\partial x_j} $$
$ u_i=\bar{u}_i+u'_i $を代入すると
$$\frac{\partial (\bar{u}_i+u'_i)}{\partial t}+(\overline{u}_j+u'_j)\frac{\partial (\bar{u}_i+u'_i)}{\partial x_j}=-\frac{1}{\rho}\frac{\partial p}{\partial x_i}+\nu \frac{\partial^2 (\bar{u}_i+u'_i)}{\partial x_j\partial x_j} $$
レイノルズ平均を取る．
$$\frac{\partial \bar{u}'_i}{\partial t}=0$$

随時更新

テンソルを用いたナビエストークス方程式の導出

テンソルを使ってナビエストークス方程式の導出手順についてメモします．

①
応力テンソル＝圧力＋粘性応力テンソル
応力テンソルの対角成分は面に垂直にかかる力なので応力テンソルのトレースの平均を圧力とする(外向き正)
つまり，応力テンソル＝平均応力テンソル＋偏差応力テンソルを考えている

$$
\mbox{応力テンソル　}\sigma= \left(
    \begin{array}{ccc}
      -p & 0 & 0 \\
      0 & -p & 0 \\
      0& 0 & -p
    \end{array}
  \right) +\mbox{粘性応力テンソル　}\tau
$$

②
ニュートン流体は偏差応力テンソル（粘性応力テンソル，せん断応力)が変形速度テンソル(ひずみ率テンソル)の線形結合で表せる　（剛体回転するときの応力は圧力のみ）

速度勾配テンソル${\bm{l}}$
$$ l_{ij}=\frac{\partial u_i}{\partial x_j} $$
直交分解する
$ {\bm{d}} $を変形速度テンソル(ひずみ速度)，${\bm{w}}$を回転速度テンソル（スピンテンソル）として$$ {\bm{d}}=\frac{{\bm{l}}+{\bm{l}}^t}{2} $$
成分は
$$ d_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right) $$

となり$$ {\bm{w}}=\frac{{\bm{l}}-{\bm{l}}^t}{2} $$である
$ {\bm{d}} $は引き延ばしを，${\bm{w}}$ は回転を表す

流体は座標変換しても変わらないという等方性を仮定すると，変形速度テンソルの係数は4階の等方テンソルで表せる
4階の等方テンソル$ A $は次の形でかける
$$ A_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu_1\delta_{ik}\delta_{jl}+\mu_2\delta_{il}\delta_{jk} $$
変形速度テンソルは対称なので
$$\tau_{ij}=A_{ijkl}d_{kl}$$
$$=\lambda \delta_{ij}d_{kk}+(\mu_1+\mu_2) d_{ij}$$
$$=\lambda \delta_{ij}d_{kk}+2\mu d_{ij}$$
\begin{equation}=\lambda \delta_{ij}(d_{11}+d_{22}+d_{33})+2\mu d_{ij}\end{equation}
とかける
偏差応力テンソルのトレース＝０なので$$ \tau_{11}+\tau_{22}+\tau_{33}=0$$
から(1)をつかって
$$(3\lambda+2\mu)\mbox{div}{\bm{u}}=0  $$
よってストークスの仮定を使うときは$$\lambda=-\frac{2}{3}\mu$$とおく
以上より応力テンソルは
$$\sigma_{ij}=-p\delta_{ij}-\frac{2}{3}\mu\delta_{ij}\mbox{div}{\bm{u}}+\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$
非圧縮なら連続の式から$ div u=0 $なので
$$\sigma_{ij}=-p\delta_{ij}+\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$となる

コーシーの運動方程式から(もしくはテンソルの力の向きを考えて)x成分を考えると
$$\rho \frac{Du_x}{Dt}=f_{x}+\frac{\partial \sigma_{xx}}{\partial x}+\frac{\partial \sigma_{xy}}{\partial y}+\frac{\partial \sigma_{xz}}{\partial z}$$
$$=f_{x}+\frac{\partial}{\partial x}(-p+2\mu\frac{\partial u_x}{\partial x}-\frac{2}{3}\mu\mbox{div}{\bm{u}})+\mu\frac{\partial}{\partial y}(\frac{\partial u_y}{\partial x}+\frac{\partial u_x}{\partial y})+\mu\frac{\partial}{\partial z}(\frac{\partial u_z}{\partial x}+\frac{\partial u_x}{\partial z})$$

$$ =f_{x}-\frac{\partial p}{\partial x}+\mu(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2})u_x+\frac{\mu}{3}\frac{\partial}{\partial x}\mbox{div}{\bm{u}}$$

よって導けた
非圧縮の場合はdiv u=0なので
$$\rho \frac{D{\bm{u}}}{Dt}={\bm{f}}-\nabla p+\mu \nabla^2 {\bm{u}} $$

圧縮性流体の運動方程式(運動量保存則)

（一般の）圧縮性流体の運動方程式は、テンソルが出てきて分かりにくいので、レイノルズの輸送定理を用いた具体的な導出をメモします。

レイノルズの輸送定理
$$ \int_{B_t}f({\bm{x}},t) dv $$
空間座標でtの変化とともに積分領域が変化する時に上の物質微分をしたいとすると、積分領域が変化するのでそのままでは上手くいかない。
一度物質座標に変換（変形勾配テンソルを使う）すれば、物質座標では積分領域が変化しない。なぜなら、tに合わせて空間座標における物質座標も変化するからである。（ゴムをのばすとゴム上の座標ものびるイメージ）
なので、物質座標で積分してから、また空間座標に戻せば良い。すると以下の結果が得られる。
$$
\frac{D}{Dt}\int_{B_t}f({\bm{x},t}) dv= \int_{B_t}(\dot{f}({\bm{x},t})+f({\bm{x},t})\mbox{div}{\bm{v}})dv
$$
$ f $の物質微分を考えて、
$$
\frac{D}{Dt}\int_{B_t}f({\bm{x},t}) dv= \int_{B_t}(\frac{\partial f}{\partial t}+{\bm{v}}\cdot\mbox{grad}f+f({\bm{x},t})\mbox{div}{\bm{v}})dv $$
よって
$$ \frac{D}{Dt}\int_{B_t}f({\bm{x},t}) dv=\int_{B_t}(\frac{\partial f}{\partial t}+\mbox{div}(f{\bm{v}})) dv
$$

運動量変化は以下のように表される
$$ \frac{D}{Dt}\int_{V}\rho{\bm{v}}dv $$
速さがベクトルだと扱いづらいので$ {\bm{e_k}} $成分について考える。レイノルズの輸送定理をつかうと
\begin{equation} \frac{D}{Dt}\int_{V}\rho{v_k}dv=\int_V(\frac{\partial \rho v_k}{\partial t}+\frac{\partial \rho v_k v_j}{\partial x_j})dv\end{equation}

また、$$ \nabla\cdot(\rho {\bm{v}}\otimes{\bm{v}})= ({\bm{e_i}}\frac{\partial }{\partial x_i})\cdot (\rho v_j v_k{\bm{e_j}}\otimes{\bm{e_k}})$$
$$=\frac{\partial \rho v_j v_k}{\partial x_j}{\bm{e}}_k$$
は(1)の右辺第二項と一致する
よって圧縮性流体の運動方程式の右辺は
$$
\frac{\partial\rho{\bm{u}}}{\partial t}+\nabla\cdot(\rho{\bm{v}}\otimes{\bm{v}})=\mbox{外部からの力積}
$$