与えられた式は鵜呑みにしないほうがいい？

前回の記事から一ヶ月。もう月日が進むのが速いね・・・。

前回は、プリミティブ方程式の話をしていたわけだが、その後に気づく。それだけではいけなかったことを。

ところで、車が走っているところを想像してみてください。この車がある速度で走っている時、この車がある時間でどこかまで走ったとしましょう。

小中では、「速度×時間＝道のり」などと習うわけですが、現実世界では、それがそのまま通用しないことは、なんとなくわかってもらえるのではないかと思います。なぜなら、他の車がいたら、場合によっては渋滞にハマるかもしれない。そうでなくても、信号で赤になっちゃってると、止まらないといけませんね。

時速60kmで1時間走れば、計算上は60km進むことになりますが、上記のようにそうならないことが多いはずです。だから、間に合いそうもないから猛スピードで飛ばせば速く着く、ってのも、実際そんなこともそれほど多く当てはまらないように思います。じゃあ、この式は役にたたないのか？もちろんそんなことはないはずです。では、どうしたら、できるだけ正確に進む距離を知ることができるだろうか、というのが今回の備忘録。説明は簡単のために、1次元で考えてみますが、実際はそうではありません。言うまでもないですが。

では、1時間を6等分して、10分ずつに分けて考えましょう。そして、そのたびに、その時の状況から改めて考え直せば、もう少し実際の所要時間に近くなるかもしれません。しかし、交通状況などでずれてしまう可能性もあります。ですが、この方法なら、最初の方法より断然確度は高くなるだろうことは予想できるのではないでしょうか。さらに言えば、その時間間隔の取り方も短いともっと確度が上がるかもしれませんね。

コンピュータで同じようなことをやってもらおうと思ったらどうでしょうか。速さ・時間・道のりの式と、1時間を小分けにして計算させたところで、うまく予想できるでしょうか。たとえば、それぞれの10分のところで速さが変わったことをプログラムはどうやって知ればいいでしょうか。人間は、暗算なり何なりでその都度臨機応変に計算して求めることができますが、コンピュータには、なんかしらの形でその情報を与えないとうまく予想してくれません。そこで、その考え方の1つにアジョイント法(もしくは四次元変分法)というのがあります。

自分は数学の専門ではなかったので、数学的な細かい話はできませんが、この方法は、計算の予想値と実測値を組み合わせて、その情報から、信頼度の高そうなデータ値を推定することでできるだけありうる値を求めるという、そういう手法です。ちなみに、上記の話は少しでもわかりやすくなるように、多少端折ったりしてるところがあります。

長くなってしまったので、ちょっとだけ実践的？な話は次回に続く。

ハエって速ぇ～速度のこととかナブラとか

速度って、一般的には、ある大きさを持ち、ある方向を向いているベクトル量だと言われる。

ところで、先の記事で書いたように、プリミティブ方程式について調べていると、調べることが山のように。自分にとってエベレスト級のつもり方。

その中で、興味深かったのが、速度の話。速度には二種類あるんだそうだ。無論、絶対速度と相対速度ではないよ。

普通、車にはスピードメータがあって、速度がわかるようになっている。ところで、風速って話はどうだろうか。風は流れていることはわかってもらえると思うが、でも、どうやって調べているのか？そりゃあんた、風速計使ってるに決まってるじゃないか、と言われればそうなのだが、よくよく考えてみると、車の場合って、車の位置がある時間 $\Delta t$ の間に動いたその変化量、いわば「位置変化量」を考えているのに対して、風速を調べる風見鶏とかって、屋根に据え付けられていて、風見鶏は動かない。位置が動かないのに位置変化量？それはないんじゃないの、と思うわけだが、それでも、実際には風速って観測されている。これらを流体力学の世界では、前者を「ラグランジュ速度」といい、後者を「オイラー速度」と呼ぶらしい。自分は物理系を専攻していなかったので、流体力学ははっきり言って何言ってのかさっぱり。記号がもうね・・・。でも grad (勾配:gradient)だけは知らないとプログラムが作れなさそうな気がしたので、ちょっとだけ必死になった。でも、grad の前に $\nabla$ を調べないと・・・。
ということで、
$\nabla$ を「ナブラ」というようです。これ、ベクトルらしいです。
$\nabla = (\cfrac{\partial}{\partial x_1}, \cfrac{\partial}{\partial x_2}, \cfrac{\partial}{\partial x_3})$ ってことらしい。単純な自分は分子(正確には分子ではないはずだ)が $\partial$ の形になっていて、「ん、何を偏微分すりゃいいんだ？」と悩みまくった(マジ)。スカラーを掛ければいいんだ、と言われれば、ベクトルのスカラー倍のことくらいは知っているが。「スカラー関数をかける」とか言われるともはやパニック。
$\nabla = \vec{i}\cfrac{\partial}{\partial x} + \vec{j}\cfrac{\partial}{\partial y} + \vec{k}\cfrac{\partial}{\partial z}(ただし、\vec{i}, \vec{j}, \vec{k}は x, y, z 方向の単位ベクトル)$ とか言われて、何がなんだか・・・。
※自分は、ベクトルをボールド？で書くとかってわかりづらくて学生のときから嫌いだったので、ベクトルの頭に矢印つけた表し方にしています。じゃあ、ナブラにも矢印つけれ、ってのは言わないでね。

しかし、まだこれ以外にもいろいろと苦戦したところがあったので、格闘の末に得たものはまた別の記事にして、後で調べやすいようにしておくことにしよう。

もう７月。

前の記事を書いたあと、プリミティブ方程式を解くプログラムに苦戦しまくって、気がついたら日付変更線またいでしまった。
ちなみに、プリミティブ方程式とは、
1) $\cfrac{\partial u}{\partial t} = -u \cfrac{\partial u}{\partial x} -v \cfrac{\partial u}{\partial y} -w \cfrac{\partial u}{\partial z} + 2 \Omega sin \phi v - \cfrac{1}{\rho}\cdot\cfrac{\partial p}{\partial x} + F_x$
2) $\cfrac{\partial v}{\partial t} = -u \cfrac{\partial v}{\partial x} -v \cfrac{\partial v}{\partial y} -w \cfrac{\partial v}{\partial z} + 2 \Omega sin \phi u - \cfrac{1}{\rho}\cdot\cfrac{\partial p}{\partial y} + F_y$
3) $g = -\cfrac{1}{\rho}\cdot\cfrac{\partial p}{\partial z}$
4) $\cfrac{\partial \rho}{\partial t} = -u \cfrac{\partial \rho}{\partial x} -v \cfrac{\partial \rho}{\partial y} -w \cfrac{\partial \rho}{\partial z} - \rho(\cfrac{\partial u}{\partial x} + \cfrac{\partial v}{\partial y} + \cfrac{\partial w}{\partial z})$
5) $\cfrac{\partial \theta}{\partial t} = -u \cfrac{\partial \theta}{\partial x} -v \cfrac{\partial \theta}{\partial y} -w \cfrac{\partial \theta}{\partial z} + H$
6) $\cfrac{\partial q}{\partial t} = -u \cfrac{\partial q}{\partial x} -v \cfrac{\partial q}{\partial y} -w \cfrac{\partial q}{\partial z} + M $
7) $P = \rho RT$
という７つの式のことを指しています。気象予報のためのプログラムのおおもとということになるのでしょう。
1), 2) は東西方向と南北方向のそれぞれの風の運動方程式です。形から、風速の時間変化率を算出しています。
3) は鉛直方向の風の運動方程式です。一般的には、およそ静力学平衡の式が成り立つと考えて計算されることが多いですが、局地的な現象は静力学平衡の式では求められませんので、非静力学平衡の式で計算される、またはオメガ方程式と呼ばれる方程式で求められるようです。自分は、3) の$g$ (重力加速度)を移行して、0 を返すようにしてみました。
4) は連続の式、もしくは質量保存の式を指します。ある体積の場所に吹き込む空気と吹き出す空気の量が釣り合うというものだったと記憶しています。
5) は熱力学方程式、もしくは熱エネルギー保存の式を指します。ここで、$\theta$ は温位と呼ばれる、空気の断熱変化が起こっても保存される、温度と似て非なるものを指しています。最後の $H$ は非断熱過程による加熱を指す量ということです。
6) は水蒸気の輸送方程式または水蒸気保存の式を指します。空気と同様水蒸気でも同じことが言えるということなのでしょう。$ M $ は非断熱過程による加湿を指す量になります。
7) は高校でも習った気体の状態方程式です。高校では $PV = (n)RT$ で習っているかもしれませんが、空気の分子量が 28.96 だとすると、これで n を割って質量が計算できます。($\cfrac{n}{28.96} = m (m: 質量)$) また、密度は質量を体積で割ったもの ($\rho = m / V$)ですから、結果として 7) 式を得ることができます。この式では、温度(温位)と密度から圧力を求めるという計算をしています。
Web 見てて、簡単なプログラムがあったので、それに触発されて作ったのですが、どうせ作るなら、できるだけそれっぽいのを、と思っていたらこんなに時間がかかってしまった・・・。

一般気象学

しかし、ここに来るまでに、いろいろとわからないことが多かったから、勉強する時間も多く使ってしまったが、ちょっとは役にたったんだから、これでよしとしよう。
今後も、自分が学んだりしたことをメモしていきます。・・・あくまでも自分のための備忘録ですが・・・。

たとえば複素数

高校で出てくる「あい」がある奴です。
数字の世界は、拡張していくと、いちばん大きい集合が複素数らしい。
ところで、複素数ってのは、虚数単位 $i$ (物理では $i$ は電流をさすので、$j$ を使うらしい)が出てくるが、それをどうやって扱おうか、という話。
ちなみに、ある複素数 $x$ において、実部 ($a + bi$ の $a$ のほう)は、共役複素数 $\bar{x}$ を使って Re$(x) = \cfrac{x + \bar{x}}{2}$ と、虚部は Im$(x) = \cfrac{x - \bar{x}}{2}$ とすれば求めることができるが、それを、プログラムのことあるごとに使っていたらかなり面倒だ。そこで、複素数を構造体にしてしまえば、扱いが楽になりそうだ。もっとも、人が扱う数学では、分数が一般でも、プログラムではそうもいかないので、とりあえず、小数で扱うことにする。そうすれば、あとは簡単だろう。

struct complex
{
    double real;
    double imaginary;
}
typedef struct complex complex;

さらに、typedef しておくと便利だ。
なお、複素数の積は、原点と複素数が表す点の長さを何倍かして、ある角度だけ回転させるものを意味する。でも、一般的にはベクトルか行列でそういう処理はするのかもしれない。

小数を分数になおす

さて、前回の記事にて、
programandmath.hatenablog.com

ってのを書きましたが。今回は、小数を分数になおす方法。
ちなみに、$\cfrac{f}{i} \cdot i = f$ が成り立つことがわかると思うが、このことから、ある小数($\cfrac{f}{i}$)に、任意の $i$ を掛ければ、整数 $f$ になるということを利用する。
したがって、小数 $\cfrac{f}{i} (= n)$ を $i$ 倍したものをある整数 $f_1$ とし、これに 1 を加えたものを $f_2$ とする。
ちなみに、そんな $i$ をどう決めるのか？ これは、for ループを使えばよい。つまり、ある範囲の中から、ループを回していって、うまく $f$ が整数になるものを探すわけだ。
ところで、これで本当にその $i$ が求めれば、該当する分数になるかというと、そうならないかもしれない。そこで、この $i$ で、$f_1$ と $f_2$ を実際に $i$ で割ってみれば良い。そして、その値のどちらかが、入力値 $\cfrac{f}{i} (= n)$ と、ある程度の(無視してもよいような)誤差範囲内にあれば、結果として、それは求める分数とみなしてもよいということになる。
以下、実際のコード。いろいろと参考にしたページがあったのを、自分で C言語で組み直した感じ。

void dec2frac(double n, long max, int diff, long frac[])
{
         long i, j1, j2;
         double v1, v2;
         for(i = 1; i <= max; i++)
         {
                 j1 = (long)(n * i); j2 = j1 + 1;
                 v1 = (double)j1 / (double)i;
                 v2 = (double)j2 / (double)i;

                 if(fabs(v1 - n) < pow(10, -(diff)))
                 {
                         frac[0] = j1; frac[1] = i;
                         return;
                 }
                 else if(fabs(v2 - n) < pow(10, -(diff)))
                 {
                         frac[0] = j2; frac[1] = i;
                         return;

                 }
         }
         frac[0] = 0; frac[1] = 0;
         return;
 }

整数と小数との判別

今までプログラムを組んでて、特に使う必要性を感じなかったのだが、調べてみて納得してしまったので、扱うことにした。

一般に小数では、切り上げと切り捨てによる整数化することができると思うのだが、切り上げ、つまり ceil() と切り捨て floor() を使うと、整数なら一致するはずだが、小数の場合、切り上げしたものと切り捨てしたものの間で 1 の差が出ることになる。なので、

    if(ceil(x) == floor(x))
        printf("%f は整数です。\n", x);
    else
        printf("%f は整数ではありません。\n", x);

とすることで、整数と小数が分けられるという。
後で、小数を分数に直す処理についても書きたいと思う。

このブログについて

このブログは、自分がプログラムだったりその他諸々の勉強して得たことを、備忘録として公開するものです。

今後公開していく予定のことがらについて、参考になさる方がいらっしゃるならば、ご自由に参考にしてくださって構いません。リンクなどもご自由にどうぞ。

ですが、もともとが自分が勉強したことについて載せるもので、間違い等もあるだろうと思います。ゆえに、参考になさる際は、自己責任にてお願いいたします。