https://beta.atcoder.jp/contests/abc044/tasks/arc060_a
すべてのカードからaを引く
平均がa⇔和がa
になる
mapを使って各カードが何枚ずつ有るか記録する。
カードがn種類あるとすると、
0<=i<=nとして、
dp[i][j]:=i番目までのカードを使って和がjになる選び方の総数
とする。
dp[0][0]=1である。
i+1番目のカードの番号がkとし、m枚あるとすると
0<=l<=mとして、i+1番目のカードをl枚使うと
dp[i+1][j+l*k]+=dp[i][j]＊mCl (コンビネーション)
になるので解ける
覚えること：
数え上げでわからない→何かを固定してdp
AtCoder Grand Contest 026
最近、競プロ始めたのでメモ。
B問題でgcdを使う部分がわからなくて理由を考えたのでメモ
あってるかどうかはわかりません。
$ g=GCD(B,D) $とする。$ B=gb, D=gd $とする時、$b,d$は互いに素。すると、
$ m=0,1,\cdots,b-1$に対して、$ dm\ \mbox{mod} \ b $は$0,1,\cdots, b-1$になる。この理由は以下の通りに考えることができる。
理由
背理法を用いる。もし、$0,1,\cdots, b-1$にならないとすると、値がかぶるので、
$$ dm_1\ \mbox{mod} \ b =x $$
$$ dm_2\ \mbox{mod} \ b =x $$
$ m_1 \neq m_2 $なので、$$ d(m_1-m_2)=0 \ \mbox{mod} \ b $$
$ (m_1-m_2) $は0より大きくb未満なので、互いに素に矛盾。
よって示せる。
一度アルゴリズムとかを真面目に勉強しないときつそうです
就活までに1600超え目指したいですが、まずは9月に1200超え目指します。
B問題でgcdを使う部分がわからなくて理由を考えたのでメモ
あってるかどうかはわかりません。
$ g=GCD(B,D) $とする。$ B=gb, D=gd $とする時、$b,d$は互いに素。すると、
$ m=0,1,\cdots,b-1$に対して、$ dm\ \mbox{mod} \ b $は$0,1,\cdots, b-1$になる。この理由は以下の通りに考えることができる。
理由
背理法を用いる。もし、$0,1,\cdots, b-1$にならないとすると、値がかぶるので、
$$ dm_1\ \mbox{mod} \ b =x $$
$$ dm_2\ \mbox{mod} \ b =x $$
$ m_1 \neq m_2 $なので、$$ d(m_1-m_2)=0 \ \mbox{mod} \ b $$
$ (m_1-m_2) $は0より大きくb未満なので、互いに素に矛盾。
よって示せる。
一度アルゴリズムとかを真面目に勉強しないときつそうです
就活までに1600超え目指したいですが、まずは9月に1200超え目指します。
gitlab自分用メモ
gitlabのメモ
超初心者なので間違えている可能性あり
レポジトリをもってくる
git clone git@gitlab.com:USERNAME/NAME.git
ローカルで持ってきたレポジトリのディレクトリにうつればgitが使える
ファイルをアップロード
git add file.cpp
git commit -m "コメント"
git push
更新
クローンしたディレクトリに入り、
git pull
消す
git rm a.txt
git commit -m "delete"
git push
ローカルをリモートに合わせる
git fetch origin
git reset --hard origin/master
設定を見る
git config --list
状態を見る
git status
超初心者なので間違えている可能性あり
レポジトリをもってくる
git clone git@gitlab.com:USERNAME/NAME.git
ローカルで持ってきたレポジトリのディレクトリにうつればgitが使える
ファイルをアップロード
git add file.cpp
git commit -m "コメント"
git push
更新
クローンしたディレクトリに入り、
git pull
消す
git rm a.txt
git commit -m "delete"
git push
ローカルをリモートに合わせる
git fetch origin
git reset --hard origin/master
設定を見る
git config --list
状態を見る
git status
C++でMPIを使ってノンブロッキング通信をする
MPIを使ってみました。
自分のパソコンは複数ノードないのでopenmpで十分な気もしますが。
簡単にまとめると
MPI：プログラムを複数プロセスで並列化する。並列する際にメモリは共有されない
Openmp：プロセスを複数スレッドに分けて並列化する。メモリは共有する
例えば1ノード4コアで4ノードあるとします。全部で16コアあるので、もし8プロセスで並列化した場合は1プロセスあたり2コアが割り当てられます。もし、openmpで1プロセス2スレッドで並列化した場合は結局1スレッド1コア使えることになります。1スレッドあたりのコアが1未満だったらどうなるのかはよく知りません。
statusはMPI_SOURCE, MPI_TAG, MPI_ERRORからなる構造体で stat.Get_source(); stat.Set_source(); stat.Get_tag(); stat.Set_tag(); stat.Get_error(); stat.Set_error();で取得できるみたいですが、Irecvを使ったときはget_source()を使っても-1になってしまい、うまくいきませんでした。
0~12000を3分割して各プロセスで和を求め、rank=0のプロセスに送信するプログラムです。
#include<iostream>
#include<mpi.h>
using namespace std;
int sum(int a, int b) {
int x = 0;
for (int i = a; i < b; i++) {
x += i;
}
return x;
}
int main(int argc, char** argv) {
MPI::Init(argc, argv);
int n = 12000;
int size, num, myid, x;
int tag = 1;
size = MPI::COMM_WORLD.Get_size();
myid = MPI::COMM_WORLD.Get_rank();
MPI::Status* st = new MPI::Status[size];
MPI::Request* req = new MPI::Request[size];
num = n / size;
if (myid == 0) {
// statusのエラーを1にする
for (int i = 0; i < size - 1; i++) {
st[i].Set_error(1);
cout << "myid="<<i<<" error before " << st[i].Get_error() << endl;
}
int* buf = new int[size];
x= sum(0, num);
cout << "myid=" << myid << " x=" << x << endl;
buf[myid] = x;
// 各プロセスから結果を受信し、配列bufに格納
for (int i = 1; i < size; i++) {
req[i-1] = MPI::COMM_WORLD.Irecv(buf + i, 1, MPI::INT, i, tag);
}
// waitする前にtestallを試す
cout << "recv flag before "<<req[0].Testall(size-1,req)<<endl;
req[0].Waitall(size - 1, req, st);
cout << "recv flag after " << req[0].Testall(size - 1, req) << endl;
// statusのエラーが0になっているか確認
for (int i = 0; i < size - 1; i++) {
cout << "myid=" << i << " error after " << st[i].Get_error() << endl;
}
// bufの結果を確認
cout << "result " << endl;
for (int i = 0; i < size; i++) {
cout << "buf" << i << "=" << buf[i] << endl;
}
}
else {
if (myid != size - 1) {
x = sum(num*myid, (myid + 1)*num);
cout << "myid=" << myid << " x=" << x << endl;
req[myid - 1] = MPI::COMM_WORLD.Isend(&x,1, MPI::INT, 0, tag);
}
else {
x = sum(num*(size-1), n+1);
cout << "myid=" << myid << " x=" << x << endl;
req[myid - 1] = MPI::COMM_WORLD.Isend(&x, 1, MPI::INT, 0, tag);
}
req[0].Waitall(size - 1, req, st);
}
MPI::Finalize();
return 0;
}
実行結果
mpic++ ex.cpp
mpirun -np 3 ./a.out
myid=0 error before 1
myid=1 error before 1
myid=0 x=7998000
myid=1 x=23998000
myid=2 x=40010000
recv flag before 0
recv flag after 1
myid=0 error after 0
myid=1 error after 0
result
buf0=7998000
buf1=23998000
buf2=40010000
自分のパソコンは複数ノードないのでopenmpで十分な気もしますが。
簡単にまとめると
MPI：プログラムを複数プロセスで並列化する。並列する際にメモリは共有されない
Openmp：プロセスを複数スレッドに分けて並列化する。メモリは共有する
例えば1ノード4コアで4ノードあるとします。全部で16コアあるので、もし8プロセスで並列化した場合は1プロセスあたり2コアが割り当てられます。もし、openmpで1プロセス2スレッドで並列化した場合は結局1スレッド1コア使えることになります。1スレッドあたりのコアが1未満だったらどうなるのかはよく知りません。
statusはMPI_SOURCE, MPI_TAG, MPI_ERRORからなる構造体で stat.Get_source(); stat.Set_source(); stat.Get_tag(); stat.Set_tag(); stat.Get_error(); stat.Set_error();で取得できるみたいですが、Irecvを使ったときはget_source()を使っても-1になってしまい、うまくいきませんでした。
0~12000を3分割して各プロセスで和を求め、rank=0のプロセスに送信するプログラムです。
#include<iostream>
#include<mpi.h>
using namespace std;
int sum(int a, int b) {
int x = 0;
for (int i = a; i < b; i++) {
x += i;
}
return x;
}
int main(int argc, char** argv) {
MPI::Init(argc, argv);
int n = 12000;
int size, num, myid, x;
int tag = 1;
size = MPI::COMM_WORLD.Get_size();
myid = MPI::COMM_WORLD.Get_rank();
MPI::Status* st = new MPI::Status[size];
MPI::Request* req = new MPI::Request[size];
num = n / size;
if (myid == 0) {
// statusのエラーを1にする
for (int i = 0; i < size - 1; i++) {
st[i].Set_error(1);
cout << "myid="<<i<<" error before " << st[i].Get_error() << endl;
}
int* buf = new int[size];
x= sum(0, num);
cout << "myid=" << myid << " x=" << x << endl;
buf[myid] = x;
// 各プロセスから結果を受信し、配列bufに格納
for (int i = 1; i < size; i++) {
req[i-1] = MPI::COMM_WORLD.Irecv(buf + i, 1, MPI::INT, i, tag);
}
// waitする前にtestallを試す
cout << "recv flag before "<<req[0].Testall(size-1,req)<<endl;
req[0].Waitall(size - 1, req, st);
cout << "recv flag after " << req[0].Testall(size - 1, req) << endl;
// statusのエラーが0になっているか確認
for (int i = 0; i < size - 1; i++) {
cout << "myid=" << i << " error after " << st[i].Get_error() << endl;
}
// bufの結果を確認
cout << "result " << endl;
for (int i = 0; i < size; i++) {
cout << "buf" << i << "=" << buf[i] << endl;
}
}
else {
if (myid != size - 1) {
x = sum(num*myid, (myid + 1)*num);
cout << "myid=" << myid << " x=" << x << endl;
req[myid - 1] = MPI::COMM_WORLD.Isend(&x,1, MPI::INT, 0, tag);
}
else {
x = sum(num*(size-1), n+1);
cout << "myid=" << myid << " x=" << x << endl;
req[myid - 1] = MPI::COMM_WORLD.Isend(&x, 1, MPI::INT, 0, tag);
}
req[0].Waitall(size - 1, req, st);
}
MPI::Finalize();
return 0;
}
mpic++ ex.cpp
mpirun -np 3 ./a.out
myid=0 error before 1
myid=1 error before 1
myid=0 x=7998000
myid=1 x=23998000
myid=2 x=40010000
recv flag before 0
recv flag after 1
myid=0 error after 0
myid=1 error after 0
result
buf0=7998000
buf1=23998000
buf2=40010000
有限要素法(FEM)の概略をメモ
重み付き残差法の大まかな流れをメモする。
ここでは
$$ \frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f=0 $$を有限要素法で解くことを考える
これは、任意の重み関数$ w $をかけて領域積分した値が0に等しいという以下の式に等価である。
$$ \int_{\Omega}w\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f \right)dS=0 $$
このまま($w$が任意)では難しいので$w=w_i$ $i=1,2\cdots ,n $について満たす解を求める
領域を少領域に分け、それぞれの領域で（例えば）1次式で解を
$$ u\approx \alpha_1+\alpha_2x+\alpha_3y $$と近似する。ただし、領域ごとに係数は異なる。このままでは扱いにくいので領域に接する節点における物理量で補間して解を表すことを考える。
小領域が三角形とすると、小領域$ e $、接点$j=(x_j,y_j)$における物理量は
$$ u^e_j=\alpha_1+\alpha_2x^e_j+\alpha_3y^e_j $$
なので、
$$
\begin{pmatrix}
u^e_1 \\
u^e_2 \\
u^e_3 \\
\end{pmatrix}
=
\begin{pmatrix}
1&x^e_1&y^e_1 \\
1&x^e_2&y^e_2 \\
1&x^e_3&y^e_3 \\
\end{pmatrix}
\begin{pmatrix}
\alpha_1\\
\alpha_2\\
\alpha_3\\
\end{pmatrix}
$$
となる。
これを解くと、$\alpha_i$が求めることができ、解（の近似）が求められる。つまり、$u^e_j$が求まれば、解が求まることがわかる。
扱いやすいように形を整えるために、
$$N^e_j(x,y)=a^e_j+b^e_jx+c^e_jy $$を用いて以下のように表すことが出来る。ただし、$a^e_j,b^e_j,c^e_j$は$x^e_i,y^e_i$を用いて表され、既知である。
$$ u\approx=\alpha_1+\alpha_2x+\alpha_3y=\sum_{j=1}^3N^e_j(x,y)u^e_j$$
$N^e_j(x,y)$は、内挿関数や形状関数と呼ばれ、節点の物理量で補間するための重み（割合）であると考えることが出来る。
接点における重みを考えると、
$$ N^e_j(x_i,y_i)=\delta_{ij} $$
であり、任意の点$(x,y) $における解において各接点の割合の和は１になることから、
$$\sum_{j=1}^3N^e_j(x,y)=1$$
となる。
ここで一度原点に戻ると、少領域$e$で
$$ \int_{\Omega_e}w_i\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f \right)dS=0 $$
となる$u$を求めたかったわけだが、$u=\sum_{j=1}^3N^e_j(x,y)u^e_j$を代入すると
$$ \int_{\Omega_e}w_i\left(\frac{\partial^2 \sum_{j=1}^3N^e_j(x,y)u^e_j}{\partial x^2}+\frac{\partial^2 \sum_{j=1}^3N^e_j(x,y)u^e_j}{\partial y^2}+f \right)dS=0 $$
これは以下のように変形できる。
$$
\sum_{j=1}^3 \int_{\Omega_e}w_i\left(\frac{\partial^2 N^e_j}{\partial x^2}+\frac{\partial^2 N^e_j}{\partial y^2} \right)dSu^e_j+\int_{\Omega_e}w_ifdS=0
$$
ここからガラーキン法(Galerkin法)を使う。
$$ w_i=-N^e_i(x,y) $$とする。ただし、小領域外では0であり、マイナスは式変形しやすくするために便宜上つけたものである。
すると、
$$
\sum_{j=1}^3 \int_{\Omega_e}\left(-N^e_i(x,y)\frac{\partial^2 N^e_j}{\partial x^2}+-N^e_i(x,y)\frac{\partial^2 N^e_j}{\partial y^2} \right)dSu^e_j=\int_{\Omega_e}N^e_ifdS
$$
これを式1とする。
ここで、次に行う式変形のために2次元のガウスの発散定理を思い出す。
$$ \int \mbox{div}{\bf{F}} dS=\int {\bf{F}}\cdot {\bf{n}}dC $$
$ {\bf{F}}=\left(N_i\frac{\partial N_j}{\partial x},N_i\frac{\partial N_j}{\partial y}\right)$を代入すると
$$\int\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+N_i\frac{\partial^2 N_j}{\partial x^2}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}+N_i\frac{\partial^2 N_j}{\partial y^2}dS=\int\left(N_i\frac{\partial N_j}{\partial x}n_x+N_i\frac{\partial N_j}{\partial y}n_y\right)dC$$
これを用いると式1は（面倒なので$e$省略）
$$
\sum_{j=1}^3 \int_{\Omega_e}\left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}\right)dSu^e_j=\int_{\Omega_e}N_ifdS+\int_{\partial\Omega_e}N_i\frac{\partial u}{\partial {\bf{n}}}dC
$$
ただし最後の項は法線方向微分である。
$$ k^e_{ij}=\int_{\Omega_e}\left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}\right)dS$$とすると、$N^e_j(x,y)=a^e_j+b^e_jx+c^e_jy $より、
$$ k^e_{ij}=\int_{\Omega_e}\left(b_ib_j+c_ic_j\right)dS$$なので計算できる。
$$ \frac{\partial u}{\partial {\bf{n}}} $$は境界条件を設定するのに使える。
$K_e=(k^e_{ij})$とすると、$K_e$は対称行列であり要素剛性行列と呼ばれる。
各要素の$K_e$を上手く合わせれば1次連立方程式に帰着され、$ u^e_j $を求めることが出来る。
ここでは
$$ \frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f=0 $$を有限要素法で解くことを考える
これは、任意の重み関数$ w $をかけて領域積分した値が0に等しいという以下の式に等価である。
$$ \int_{\Omega}w\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f \right)dS=0 $$
このまま($w$が任意)では難しいので$w=w_i$ $i=1,2\cdots ,n $について満たす解を求める
領域を少領域に分け、それぞれの領域で（例えば）1次式で解を
$$ u\approx \alpha_1+\alpha_2x+\alpha_3y $$と近似する。ただし、領域ごとに係数は異なる。このままでは扱いにくいので領域に接する節点における物理量で補間して解を表すことを考える。
小領域が三角形とすると、小領域$ e $、接点$j=(x_j,y_j)$における物理量は
$$ u^e_j=\alpha_1+\alpha_2x^e_j+\alpha_3y^e_j $$
なので、
$$
\begin{pmatrix}
u^e_1 \\
u^e_2 \\
u^e_3 \\
\end{pmatrix}
=
\begin{pmatrix}
1&x^e_1&y^e_1 \\
1&x^e_2&y^e_2 \\
1&x^e_3&y^e_3 \\
\end{pmatrix}
\begin{pmatrix}
\alpha_1\\
\alpha_2\\
\alpha_3\\
\end{pmatrix}
$$
となる。
これを解くと、$\alpha_i$が求めることができ、解（の近似）が求められる。つまり、$u^e_j$が求まれば、解が求まることがわかる。
扱いやすいように形を整えるために、
$$N^e_j(x,y)=a^e_j+b^e_jx+c^e_jy $$を用いて以下のように表すことが出来る。ただし、$a^e_j,b^e_j,c^e_j$は$x^e_i,y^e_i$を用いて表され、既知である。
$$ u\approx=\alpha_1+\alpha_2x+\alpha_3y=\sum_{j=1}^3N^e_j(x,y)u^e_j$$
$N^e_j(x,y)$は、内挿関数や形状関数と呼ばれ、節点の物理量で補間するための重み（割合）であると考えることが出来る。
接点における重みを考えると、
$$ N^e_j(x_i,y_i)=\delta_{ij} $$
であり、任意の点$(x,y) $における解において各接点の割合の和は１になることから、
$$\sum_{j=1}^3N^e_j(x,y)=1$$
となる。
ここで一度原点に戻ると、少領域$e$で
$$ \int_{\Omega_e}w_i\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f \right)dS=0 $$
となる$u$を求めたかったわけだが、$u=\sum_{j=1}^3N^e_j(x,y)u^e_j$を代入すると
$$ \int_{\Omega_e}w_i\left(\frac{\partial^2 \sum_{j=1}^3N^e_j(x,y)u^e_j}{\partial x^2}+\frac{\partial^2 \sum_{j=1}^3N^e_j(x,y)u^e_j}{\partial y^2}+f \right)dS=0 $$
これは以下のように変形できる。
$$
\sum_{j=1}^3 \int_{\Omega_e}w_i\left(\frac{\partial^2 N^e_j}{\partial x^2}+\frac{\partial^2 N^e_j}{\partial y^2} \right)dSu^e_j+\int_{\Omega_e}w_ifdS=0
$$
ここからガラーキン法(Galerkin法)を使う。
$$ w_i=-N^e_i(x,y) $$とする。ただし、小領域外では0であり、マイナスは式変形しやすくするために便宜上つけたものである。
すると、
$$
\sum_{j=1}^3 \int_{\Omega_e}\left(-N^e_i(x,y)\frac{\partial^2 N^e_j}{\partial x^2}+-N^e_i(x,y)\frac{\partial^2 N^e_j}{\partial y^2} \right)dSu^e_j=\int_{\Omega_e}N^e_ifdS
$$
これを式1とする。
ここで、次に行う式変形のために2次元のガウスの発散定理を思い出す。
$$ \int \mbox{div}{\bf{F}} dS=\int {\bf{F}}\cdot {\bf{n}}dC $$
$ {\bf{F}}=\left(N_i\frac{\partial N_j}{\partial x},N_i\frac{\partial N_j}{\partial y}\right)$を代入すると
$$\int\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+N_i\frac{\partial^2 N_j}{\partial x^2}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}+N_i\frac{\partial^2 N_j}{\partial y^2}dS=\int\left(N_i\frac{\partial N_j}{\partial x}n_x+N_i\frac{\partial N_j}{\partial y}n_y\right)dC$$
これを用いると式1は（面倒なので$e$省略）
$$
\sum_{j=1}^3 \int_{\Omega_e}\left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}\right)dSu^e_j=\int_{\Omega_e}N_ifdS+\int_{\partial\Omega_e}N_i\frac{\partial u}{\partial {\bf{n}}}dC
$$
ただし最後の項は法線方向微分である。
$$ k^e_{ij}=\int_{\Omega_e}\left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}\right)dS$$とすると、$N^e_j(x,y)=a^e_j+b^e_jx+c^e_jy $より、
$$ k^e_{ij}=\int_{\Omega_e}\left(b_ib_j+c_ic_j\right)dS$$なので計算できる。
$$ \frac{\partial u}{\partial {\bf{n}}} $$は境界条件を設定するのに使える。
$K_e=(k^e_{ij})$とすると、$K_e$は対称行列であり要素剛性行列と呼ばれる。
各要素の$K_e$を上手く合わせれば1次連立方程式に帰着され、$ u^e_j $を求めることが出来る。
C++でのバイナリファイルの読み書き 自分用メモ
バイナリファイルをC++で用いる際に注意することをメモする。ただし、自分用メモであり、正しさは保証しない。
注意点1:
テキストファイルのように＞＞や＜＜は使えない（っぽい）のでwriteやreadを使う必要がある
注意点2：
writeやreadは char型のポインタとデータサイズしか引数として受け取らない（っぽい）
なので、intやdoubleを読み書きする場合は
reinterpret_cast<const char*>
を用いてint やdoubleのポインタをcharのポインタに変換する
例
#include<iostream>
#include<fstream>
using namespace std;
int main(int argc, char* argv[]) {
ofstream ofs;
ofs.open("example.bin",ios::binary);
for (int i = 0; i < 3;i++) {
ofs.write(reinterpret_cast<char *> (＆i), sizeof(i));
}
ofs.close();
ifstream ifs;
int z;
ifs.open("example.bin", ios::binary);
for (int i = 0; i < 3; i++) {
ifs.read(reinterpret_cast<char *>(＆z), sizeof(int));
cout << z<<endl;
}
ifs.close();
return 0;
}
実行結果
./a.out
0
1
2
Linuxのodコマンドで確認すると、
od example.bin
0000000 000000 000000 000001 000000 000002 000000
0000014
なので上手くいっているっぽい
注意点1:
テキストファイルのように＞＞や＜＜は使えない（っぽい）のでwriteやreadを使う必要がある
注意点2：
writeやreadは char型のポインタとデータサイズしか引数として受け取らない（っぽい）
なので、intやdoubleを読み書きする場合は
reinterpret_cast<const char*>
を用いてint やdoubleのポインタをcharのポインタに変換する
例
#include<iostream>
#include<fstream>
using namespace std;
int main(int argc, char* argv[]) {
ofstream ofs;
ofs.open("example.bin",ios::binary);
for (int i = 0; i < 3;i++) {
ofs.write(reinterpret_cast<char *> (＆i), sizeof(i));
}
ofs.close();
ifstream ifs;
int z;
ifs.open("example.bin", ios::binary);
for (int i = 0; i < 3; i++) {
ifs.read(reinterpret_cast<char *>(＆z), sizeof(int));
cout << z<<endl;
}
ifs.close();
return 0;
}
実行結果
./a.out
0
1
2
Linuxのodコマンドで確認すると、
od example.bin
0000000 000000 000000 000001 000000 000002 000000
0000014
なので上手くいっているっぽい
風上差分と数値拡散
風上差分については感覚的に流れの上流の情報を使うと説明されることが多いので数学的な説明をメモする
まず、拡散方程式について軽くメモする。
拡散方程式は、以下の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条件を満たす
・濃度が高い方から低い方へ流れる
・流れの速さは濃度の勾配に比例する
$ 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)式を用いて同様のことをした場合(風上差分ではないほうの差分)右辺が負の拡散項になるので濃度勾配に対して物理量が更に濃度勾配が大きくなる作用を持つ。そのため、安定性が悪い。
登録:
投稿 (Atom)
-
（一般の）圧縮性流体の運動方程式は、テンソルが出てきて分かりにくいので、レイノルズの輸送定理を用いた具体的な導出をメモします。 レイノルズの輸送定理 $$ \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を求めれば良い ... -
テンソルを使ってナビエストークス方程式の導出手順についてメモします． ① 応力テンソル＝圧力＋粘性応力テンソル 応力テンソルの対角成分は面に垂直にかかる力なので応力テンソルのトレースの平均を圧力とする(外向き正) つまり，応力テンソル＝平均応力テンソル＋偏差応力テンソル...