数値計算
1次元のEuler方程式をRoeの方法で、MUSCLで解くのは出来たぽいので、そのあとのTODOを考えてみる。ついでに1次元の1変数の移流方程式もMUSCLで解くことが出来たので、これも別のものに流用できそうである。 2次元の直交座標系でEuler方程式をRoeの方法でMUSC…
それなりに前から仕事の合間に圧縮性のコードを解こうとしてるんだけど、まあMUSCLが多いですから、それ使おうとしてます。で、1次元の移流方程式を解こうとしてるんですが、そもそも衝撃波の問題を解くには質量保存と、運動量保存と、エネルギー保存の3つの…
TVD凄いね。マジビビ。移流方程式でオーバーシュートがねー。けど、若干波形が崩れ気味。なのが気になる。 やっぱ使いどころを考えた方が良いのかも。 それでもTVD使ってもオーバーシュートするけどな。 そういうことでこれからshock tubeでも作ってみます。…
Fortranでも構造体が使えるらしいです。素敵です。変数が増えすぎると氏ねるので、とても良いと思います。でもCのように構造体へのポインタとかいう使い方はおっかないのでできないです。 まあでもFortranの場合はサブルーチンにぶち込んでしまえば問題ない…
固有値問題に心が折れたのでMAXIMAをインスコします。3x3の行列なんだが、形式が複雑だと氏ねます。特に固有値固有ベクトルまではなんとかなるけど、右と左の対格化行列を出す段になると結構面倒臭い。 昔の人はエロかったなあと感心する次第です。 なんでこ…
大体俺がよく出すバグはコアダンプなので、デバッガを使ってそれを追跡する方法について。 コンパイルするときに、 cc -g -O0 hoge.cにしてコンパイル。 ついでに、 $ ulimit -c unlimitedにしてコアダンプをださせるようにする。そして、a.outを実行。で、…
gccのデバッガ初めて使った。VisualStudioのは使ったことがあって、それはGUIで操作できて、とても使いやすかったので、以前gccのデバッガ使ったときにCUIで取っつきが悪かったのでスルーしてた。けど、止むに止まれぬ理由で使ってみた。 やっぱデバッガは便…
主に自分が忘れない用。 こう書くと格好いい感じだけど、mallocしてるだけである。ヒーブ領域よりもスタック領域の方がでかいらしいので、プログラムの頭で、double Hoge[NN]とかして、そんな沢山配列切れないよって言われたらmallocとか、極めていい加減な…
訳あって大量にプロットしないといけなくなったのでpopenを使ってgnuplotにデータ吐かせます。 まあ、使い方はfopenとか、何とか、open関数とかと同じらしい。 #include <stdio.h> int main () { FILE *gp; gp = popen ( "gnuplot -persist", "w"); fprintf (gp, "p s</stdio.h>…
Mathematicaでなくてもgnuplotで十分いけるらしいことが分かりました。特にMathematicaの3次元プロットは使えそうで使えない。 gnuplotでデータファイルから3次元プロットするときは、ちょっとできあがりの美しさを考えたらなんらかの形でフィッティングかけ…
gnuplotならsplot使うだけ。でもデータから曲面を出すのはgnuplotでは面倒なので、Mathematicaで。 ListPlot3Dでできる。らしい。けど、Macだと何故か上手くいかない。(x,y,z)がnあったとして、それを解釈してくれないらしい。 http://reference.wolfram.com…
pngでも出せるらしい。tiffでも出せるともっと嬉しい。 出力先をpngに指定すれば良いらしい。 gnuplot > set terminal pngにする。 ここではプロットをpng形式で出力しますよというのをやってる。画面に出力する場合は、Linuxとかならx11、Windowsならwindow…
複数の窓を開く > set t x11 3とかやると窓を4つ開けるらしい。 http://jegog.phys.nagoya-u.ac.jp/~tkonishi/cn/misc/gnuplot.html 補助メモリの表示 > set mxtics 10とするとx軸上の主メモリの間に補助メモリが10点打たれるらしい。
波動方程式の初期値問題を解くときに一番ましな方法らしい。色々とあって、誤差の増幅の度合とかを調べられるんだけど、まあそこはスルー。要するに、Fourier級数的な解があるとして、そのスペクトルの大きさを調べる的なことをやるとできる。 多分蛙跳び法…
あのね、配列全部を関数にぶん投げたかったのです。 一次元配列のときは、int A[NN];という配列があったとして、それをなんらかの関数にぶち込むときは、int hoge (int A[], int n);とかでよかった訳だが、今は訳あって3次元の配列、A[LL][MM][NN]を考えてる…
LAPACKまではaptで、lapack3とlapack3-develを落としてくれば入れられる。 CLAPACKは自分でmakeしないとダメっぽい。 http://www.view.human.nagoya-u.ac.jp/watanabe-lab/graduates/01year/hiroyuki/lab/clapack/index.html の通りにやったら動いた。 もう…
まあね、ただの直交関数展開なんですがね、覚え書き程度に。 独立変数がx, tの未知変数f(x,t)があって、それに微分作用素Lが掛かって偏微分方程式、 があるとする。 これを多項式の線形和で表すとすると、 であるとして、これを微分方程式に代入すると、 こ…
TrigExpandという関数を使うらしい。超簡単。 それにしてもなんでMathematicaさまの出す答えと俺の出すのが違うのだろうか。どこにも二次形式なんて出てこないんで場合分けなんて要らないし、簡単な代数計算のはずなのになあ。おかしいと思うので検算してみ…
conmplex.hを使うと複素数を使えます。変数の宣言は、 #include <complex.h> ... int main () { double _Complex z; int _Complex a; return 0; }とかそんな感じ。 http://www.alab.t.u-tokyo.ac.jp/~bond/doc/complex.html に詳しい。complex.hを使うときに特にオプシ</complex.h>…
大量にあるデータを一辺にFFT掛けたいっつーことで、OctaveとかMathematicaでやるのはちょっとよく分からんのでC言語でやっちまう。自分でFFTのコード書くのは愚の骨頂なので、当然パッケージを使う。っつーかさあ、波数空間でやってりゃこんな面倒なことね…
どうやらemacs上でgnuplotを動かせたりするらしい。もうここまで来るとemacsしか入ってないOSとかもありかなあとか思ってみたりする。最近はOctave動かすのもemacs上でだし、TeXもemacsでだし、これでgnuplotまでemacsでっていったら、本当にあとはメーラと…
バッドノウハウっていうのは、仕様が複雑で色々覚えないと操作もままならないシステムのことらしい。そういうことになるともうUnixなんてバッドノウハウの塊みたいなもんですね。昔の数多のLinuxディストリビューションも。PlamoとかKONDARAとか、その昔のDe…
C言語とは若干勝手が違う。なんかPythonっぽい。 例えば、数字を1から10まで画面に表示させるには、 > for n = 1:10 > disp (n) > endになる。なんかwhile文っぽい。恐らく意味としては、nが1から1つずつ増えて10になるまでとかそんな感じ。今の時点ではこん…
文字列の入力とか、連結とか、分解とか。「続きを読む」で。
適当にOctaveのコマンドを書いて、それをhoge.mというファイルに保存。そして、Octaveを起動して、 > hoge.mでできる。超簡単。また、emacsでOctaveモードを使えるようにしてあったときは、スクリプトを編集してる状態から、M-x run-octaveで動きます。 Wind…
emacsでOctaveモードを使うには、~/.emacsファイルに、 ;;; octave (autoload 'octave-mode "octave-mod" nil t) (setq auto-mode-alist (cons '("\\.m$" . octave-mode) auto-mode-alist)) (add-hook 'octave-mode-hook (lambda () (abbrev-mode 1) (auto-f…
英語で小行列はminorというらしい。 行列a(m,n)があって、1≦m≦M、1≦n≦Nのうち、1≦m≦M1、1≦n≦N1を取り出すことを考えます。但し、M1≦M、N1≦N。で、その小行列をbとすると、 > b = a (1:M1, 1:N1);詳しくは「続きを読む」で。
2つのデータaiとbiがあったとして、これの相関は、 > corrcoef (a, b)ででる。終わり。超簡単。 ここで引数を行列にして、corrcoef (A)とかすると相関係数行列というものが出てくるらしいんだけど、それの意味する所はちょっとよくわからない。漏れは相関係…
超簡単です。問題はどうやってデータをインポートするかなんだけど、普通の形式でデータが与えられていればまあ問題なかろう。 手順は、 データを読み込む fft()関数でFFT 結果をプロット ある変数hogeにデータが入ってたとしたら、それをfft(hoge)とするだ…
C言語っぽくもできるらしいんだけど、それだと遅いので、Octaveのコマンドを使った方がいいと思う。Cフラグ立てて、fscanfとかsscanfとかもいいんだけど、それよりはデータを事前にOctaveの読み込みやすい形にしといてload関数でガシャポンと読み込んだ方が…
