ラグランジュ補間による差分法の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は三角なので簡単に求められる
次はQR分解で固有値求める予定
最小二乗法のプログラムは↓
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次式で近似してみる
最小二乗法のプログラムは↓
#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が入
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$$
随時更新
速度を平均と変動成分に分けます．
レイノルズ分解
$$ 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}} $$
①
応力テンソル＝圧力＋粘性応力テンソル
応力テンソルの対角成分は面に垂直にかかる力なので応力テンソルのトレースの平均を圧力とする(外向き正)
つまり，応力テンソル＝平均応力テンソル＋偏差応力テンソルを考えている
$$
\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}} $$
登録:
投稿 (Atom)
-
乱流を扱うときに使うモデルです． 速度を平均と変動成分に分けます． レイノルズ分解 $$ u_i=\bar{u}_i+u'_i $$ 流れが定常→時間平均 一様な流れ→空間平均 一般→アンサンブル平均 アンサンブル平均 $$\bar{u}=\lim_{N...
-
（一般の）圧縮性流体の運動方程式は、テンソルが出てきて分かりにくいので、レイノルズの輸送定理を用いた具体的な導出をメモします。 レイノルズの輸送定理 $$ \int_{B_t}f({\bm{x}},t) dv $$ 空間座標でtの変化とともに積分領域が変化する時に上の物質...
-
|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を求めれば良い ...