風上差分と数値拡散

風上差分については感覚的に流れの上流の情報を使うと説明されることが多いので数学的な説明をメモする
まず、拡散方程式について軽くメモする。
拡散方程式は、以下の2条件を満たす
・濃度が高い方から低い方へ流れる
・流れの速さは濃度の勾配に比例する
$ x-\Delta x,x,x+\Delta x $における濃度$ f(x) $を考えて、$ x $に流れ込む量と流れ出す量を足し合わせると拡散方程式
$$\frac{\partial f}{\partial t}=D\frac{\partial^2 f}{\partial x^2}$$
が得られる　ただし、$D$は正のときに拡散をするが、負だと濃度勾配をさらに大きくすることになる。詳しくは負の拡散係数とかでググれば出てくる

移流方程式
$$ \frac{\partial f}{\partial t}+u\frac{\partial f}{\partial x}=0$$
を数値的に解くことを考える
テイラー展開をすると
$$ f_{i-1}=f_i-\left(\frac{\partial f}{\partial x}\right)_i\Delta x+\frac{(\Delta x)^2}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+\mathcal{O}(\Delta x^3)$$
$$ f_{i+1}=f_i+\left(\frac{\partial f}{\partial x}\right)_i\Delta x+\frac{(\Delta x)^2}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+\mathcal{O}(\Delta x^3)$$
よって
\begin{equation} u\left(\frac{\partial f}{\partial x}\right)_i=u\frac{f_i-f_{i-1}}{\Delta x}+u\frac{\Delta x}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+u\mathcal{O}(\Delta x^2)\end{equation}

\begin{equation} u\left(\frac{\partial f}{\partial x}\right)_i=u\frac{f_{i+1}-f_{i}}{\Delta x}-u\frac{\Delta x}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+u\mathcal{O}(\Delta x^2)\end{equation}
上二式の二階偏微分の項が$ u $の正負によって拡散項になったり負の拡散項になることが本質である。$u$が正の時は上の1式が風上差分となるがその時にちゃんと拡散項になっていて負の拡散項ではないことを確かめる。

1式を変形すると
$$ u\frac{f_i-f_{i-1}}{\Delta x}=u\left(\frac{\partial f}{\partial x}\right)_i-u\frac{\Delta x}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+u\mathcal{O}(\Delta x^2)  $$
移流方程式に1式を用いて移流方程式を離散化すると
$$ \frac{\partial f}{\partial t}+u\frac{f_i-f_{i-1}}{\Delta x}=0$$
を解くことになるが、これは上式を代入すると
$$ \frac{\partial f}{\partial t}+u\frac{\partial f}{\partial x}=u\frac{\Delta x}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+u\mathcal{O}(\Delta x^2)$$
これは、右辺が拡散方程式と同じで拡散項になっていることがわかる。
(2)式を用いて同様のことをした場合(風上差分ではないほうの差分)右辺が負の拡散項になるので濃度勾配に対して物理量が更に濃度勾配が大きくなる作用を持つ。そのため、安定性が悪い。

ラグランジュ補間による差分法の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

共役勾配法(CG法)の収束について

工学の様々な場面で登場する前処理付き共役勾配法。
共役勾配法は前処理をするのが重要(PCG法)とよく言われますが、その理由は省略されていることが多いのでメモしておきます。

$ A $は正定値対称行列
$$ {\bm{r}}_0={\bm{b}}-A{\bm{x}}_0 $$
$$ {\bm{p}}_{i+1}={\bm{r}}_{i+1}+\beta_i{\bm{p}}_i $$
$$ {\bm{x}}_{i+1}={\bm{x}}_i+\alpha_i{\bm{p}}_i $$
$$ {\bm{r}}_{i+1}={\bm{b}}-A{\bm{x}}_{i+1} ={\bm{b}}-A({\bm{x}}_i+\alpha_i{\bm{p}}_i )={\bm{r}}_i-\alpha_iA{\bm{p}}_i $$
$$ A{\bm{x^*}}={\bm{b}} $$
$$ {\bm{e}}_i={\bm{x^{*}}}-{\bm{x}}_i $$
$$ {\bm{e}}_{i+1}={\bm{{e}}}_i-\alpha_i{\bm{p}}_i $$
ベクトルの太字は面倒なので以下省略。
クリロフ部分空間
$$ \mathcal{K}_i=\mbox{span}\{ r_0,Ar_0,\cdots A^{i-1}r_0\}= \mbox{span}\{ Ae_0,A^2e_0,\cdots A^{i}e_0\}$$
$ i $次多項式の多項式$ P_i $を用いて(polynomial)
$$ e_i=P_i(A)e_0 $$
と書ける。アルゴリズムを書いてみて動かすとわかるが、$ P_i(0)=1 $
ここからが重要
$ A $の固有値、固有ベクトルを$ \lambda_i,u_i $とする。
$ e_0 $は以下のように表せる。$ a_i $は係数。
$$ e_0=\sum_{j=1}^n a_i u_i $$
$$ e_i=P_i(A)\sum_{j=1}^n a_j u_j =\sum_{j=1}^n a_jP_i(\lambda_j)u_j$$
$$ Ae_i=P_i(A)A\sum_{j=1}^n a_j u_j =\sum_{j=1}^n a_j\lambda_jP_i(\lambda_j)u_j$$
対称行列の固有ベクトルは直交するので
$$ (e_i,Ae_i)=\sum_{j=1}^n {a_j}^2 P_i(\lambda_j)^2\lambda_j $$
$ A $は正定値対称なので、$ (x-x^*,A(x-x^*)) $を最小にするのは$ x=x^* $のみである。
共役勾配法は各ステップで$ (e_{i+1},Ae_{i+1}) $が最小になるように$ \alpha_i $を選んでいるので、
$ \mathcal{K}_i $上で$ (e_{i},Ae_{i}) $が最小になるから、
$$ (e_i,Ae_i)=\min_{P_i}\sum_{j=1}^n {a_j}^2 P_i(\lambda_j)^2\lambda_j $$
ここで、$ 0\leq X_i,Y_i $とし、$ X_{max}=\max_{i }X_i $とすると、一般に$ x^x $
$$　\sum_i X_i Y_i\leq\sum_i X_{max} Y_i=X_{max}\sum_i Y_i $$
が成り立つ。Aは正定値なので$ X_i=P_i(\lambda_i)^2,Y_i=a_j^2\lambda_j $とすると、
\begin{eqnarray}
\sum_{j=1}^n {a_j}^2 P_i(\lambda_j)^2\lambda_j&\leq  &
\max_{\lambda_k\in{\mbox{Aの固有値の集合}}}P_i(\lambda_k)^2\sum_{j=1}^n {a_j}^2 \lambda_j\\
&\leq&\max_{\lambda\in[\lambda_{min},\lambda_{max}]}P_i(\lambda)^2\sum_{j=1}^n {a_j}^2 \lambda_j
\end{eqnarray}
よって、
\begin{eqnarray}
(e_i,Ae_i)&=&\min_{P_i}\sum_{j=1}^n {a_j}^2 P_i(\lambda)^2\lambda_j\\
&\leq&\min_{P_i}\max_{\lambda\in[\lambda_{min},\lambda_{max}]}P_i(\lambda_j)^2\sum_{j=1}^n {a_j}^2 \lambda_j\\
&=&\min_{P_i}\max_{\lambda\in[\lambda_{min},\lambda_{max}]}P_i(\lambda)^2(e_0,Ae_0)
\end{eqnarray}
上式は$ \lambda $がmin~maxまで動く時の関数の最大値を最小にするような$ P_i $を選んだ時の値を表している
最小偏差多項式を使えば収束の速さが評価できる。
$ k $次のチェビシェフ多項式を以下のように書ける
$$ T_k(x)=\cos(l\arccos x) $$
$$ \cfrac{\lambda_{min}}{\lambda_{max}}
$$
が大きいほど早いことがわかる