備忘録

Perl

2014-09-11T22:10:00.002+09:00

最近PDLを使おうかと思ってます。

いや、仕事ではないですが、せっかくperlの勉強をしたいと思っているので。

自由度が高いけど、書き手によってコードのエレガント差が変わる部分が大きい言語は結構好きかもしれません。

学習レベルによってコードが変わるのは、自分の成長度がわかっていいですね。
業務では嫌われる要因でしょうが・・・。

モダンPerlの新版はいつ出るんでしょうか。
結構楽しみにしてます。

あと、Perl入学式に参加して、
「書き方は自由です」
というのは、勇気づけられました。
もちろん、技巧をこらしても良い、ということではなく、
「汚く書いても動くし、技巧をこらしても動くけど、ちゃんと考えて書いてね」
と言われてる気がします。

それはいわば、自然言語に近いのかなと。

文章の流れを考えて、くどく「私は」「私は」と書くのは、
読みづらいし、大事な所で主語が抜けるのは、本当に意味が分からない。

話し言葉と文章では異なりますし。

そういった意味では、使い手のレベルに依存し、
高度なコードも、低度なコードも可能なので、自然言語に似ていると感じます。

さすが、ラリーさんです。

今後も発展して欲しいと思いますし、無くなることはないでしょう。
Linuxにはほとんど入ってますし。

死ぬまでにどれだけ理解出来るか、楽しみです。

カルマンフィルタ

2014-06-18T22:54:00.000+09:00

足立先生の「カルマンフィルタの基礎」を勉強させていただきます。
初心者にとっては非常に良い本だと思ってます。

その中で、例題プログラムがあり、MATLABなのですが、
ビンボーなため、持っておりません。

そこでpythonで書いてみました。
最初の例題は、スカラー変数なので、若干弄っております。
コメントはしません。
色々おこられそうなので。。。

本当に勉強になります。

import numpy as np

import pylab as pl

import random

def kf(A,B,Bu,C,Q,R,u,y,xhat,P):

xhatm = A*xhat+Bu*u

Am = A;Bm=B;Cm=C;I =1

Pm = A*P*Am + B*Q*Bm

G = Pm*C/(C*Pm*Cm+R)

xhat_new = xhatm + G*(y-C*xhatm)

P_new = (I-G*C)*Pm

return xhat_new,P_new,G

A = 1.

b = 1.

c = 1.

Q = 1.

R = 5.

N = 300

x = np.zeros(N)

y = np.zeros(N)

v = np.zeros(N)

w = np.zeros(N)

xhat = np.zeros(N)

t = np.zeros(N)

for i in range(N):

v[i] = random.normalvariate(0,Q)

w[i] = random.normalvariate(0,R)

t[i] = i

#print v[i],w[i]

y[0] = c*x[0] + w[0]

for k in range(1,N):

x[k] = A*x[k-1]+b*v[k-1]

y[k] = c*x[k] + w[k]

P = 0

xhat[0] = 0.0

for k in range(1,N):

xhat[k],P,G = kf(A,b,0,c,Q,R,0,y[k],xhat[k-1],P)

pl.plot(t,y,label='measured')

pl.plot(t,x,label='ture')

pl.plot(t,xhat,label='estimate')

pl.grid(True)

pl.legend()

pl.savefig('kf')

pl.show()

とりあえずうまくいってるようです。

Harten-Yeeをちょっと修正

2014-06-15T23:43:00.001+09:00

以前のものを多少修正しました。
もうちょっと美しく出来れば良いんですけど・・・。
1次風上も残してますが、２０行くらい少なくなってます。

#advection #################################################
# Harten Yee Upwind ########################################
############################################################
#import#####################################################
from pylab import *
def cal_fai(z):
eps = 0.0000000001
if abs(z) > eps:
return abs(z)
else:
return (z*z+eps*eps)/(2*eps)
def cal_sigma(z,dt,dx):
return 0.5*(cal_fai(z)-dt/dx*z*z)
def minmod2(x,y):
if abs(x)<=abs(y) and x*y >=0.0:
return x
elif abs(x)>abs(y) and x*y >= 0.0:
return y
elif x*y<0:
return 0.0
else:
print "Error minmode2"
# Harten Yee ##############################################
class Cell:
def __init__(self,dxl,u,f,fn,fl,fr):
self.dxl = dxl
self.u = u
self.f = f
self.fn = fn
self.fl = fl
self.fr = fr
self.df = 0.0
self.gi = 0.0
self.gamma = 0.0
self.fai = 0.0
def upwind(self,dt,cellup,celldn):
ultmp = 0.5*(u + cellup.u)
urtmp = 0.5*(celldn.u + u)
if self.u > 0:
self.fl = ultmp*cellup.f
self.fr = urtmp*self.f
else:
self.fl = ultmp*self.f
self.fr = urtmp*celldn.f
def caldf(self,celldn):
self.df = celldn.f - self.f
def calgi(self,cellup):
self.gi = minmod2(self.df,cellup.df)
def calgamma(self,dt,celldn):
if self.df != 0:
self.gamma = cal_sigma(self.u,dt,self.dxl)*(celldn.gi - self.gi)/self.df
else:
self.gamma = 0.0
def calfai(self,celldn):
self.fai = cal_sigma(self.u,dt,self.dxl)*(self.gi+celldn.gi)-cal_fai(self.u+self.gamma)*self.df
def yee(self,dt,cellup,celldn):
csdn = 0.5*(celldn.u*celldn.f+self.u*self.f)
csup = 0.5*(self.u*self.f + cellup.u*cellup.f)
cfdn = 0.5*self.fai
cfup = 0.5*cellup.fai
self.fr = csdn+cfdn
self.fl = csup+cfup
def euler(self,dt):
self.fn = self.f - dt/self.dxl*(self.fr - self.fl)
def shift(self):
self.f = self.fn
#set parameters ############################################
in_ff = open('test.txt')
imx=100
dx = 1.0
dt = 0.2
u = 0.2
print "-"*40
print "u =",u
print "dx =",dx
print "dt =",dt
print "-"*40
############################################################
#initial condition ################################
cells = []
fini = []
i = 0
for x in in_ff:
print i
a = x.split("\t")
for k in range(0,1):
print a[k]
cells = cells + [Cell(float(dx),float(u),float(a[1]),float(a[1]),0.0,0.0)]
fini = fini + [float(a[1])]
i=i+1
imx=i
print "imx=",imx
print "*"*40
print "xl,f"
for i in range(0,imx):
print i," ",cells[i].f,cells[i].u
# Time integration ################################
for n in range(0,1000):
#Harten Yee##############################
for i in range(0,imx-1):
cells[i].caldf(cells[i+1])
for i in range(2,imx-1):
cells[i].calgi(cells[i-1])
for i in range(3,imx-2):
cells[i].calgamma(dt,cells[i+1])
for i in range(2,imx-1):
cells[i].calfai(cells[i+1])
for i in range(1,imx-3):
cells[i].yee(dt,cells[i-1],cells[i+1])
for i in range(2,imx-2):
cells[i].euler(dt)
#shift#########################################
for i in range(2,imx-2):
cells[i].shift()
# Exact ###########################################
fout = [0]
fex = [0]
xl = [0]
for i in range(0,imx-1):
if i>=49 and i<=59:
fex = fex +[1.0]
else:
fex = fex +[0.0]
fout = fout + [float(cells[i].f)]
xl = xl + [float(i)]
#output############################################
ff=open("output.dat","w")
for i in range(0,imx-1):
ff.write(str(i))
ff.write(" ")
ff.write(str(fout[i]))
ff.write("\n")
ff.close()
#####################graph#########################
y1=fout
y2=fini
y3=fex
print len(fout),len(fini),len(fex)
plot(xl,y1,label = 'Harten-Yee')
plot(xl,y2,label = 'INITIAL')
plot(xl,y3,label = 'EXACT')
xlabel('x')
ylabel('f')
title('Hatren-Yee')
axis([-0,100.0,-0.2,1.2])
legend()
grid(True)
savefig('Harten-Yee')
show()
###################################################

Java

2014-05-11T16:06:00.001+09:00

今、仕事でJavaを使ってます。

正直、FortranやC的な書き方しか出来てません・・・。

今更ですが、自分自身オブジェクト指向が出来ていないと思ってます。

結城先生のJavaの本を読んで、多少わかった気にはなってますが、
まだまだだなと。

実装自体は、数値解析的な所かつライブラリを使う所なので難しい部分はありますが・・・。

Fortran

2014-03-22T09:50:00.001+09:00

新しい会社では全く使わなくなったFortran＋本の話です。

「数値計算のためのFrotran90/95プログラム入門」は以前の会社に入って良く読みました。京都大学の牛島先生の本です。

特定の分野では77は使われています。
ゼロから書く場合は90/95へのシフトがなされているようです。

大学の時は77が多くて、
・Windowsかつ小さなプログラムならfcpad＋salfordで解く
・Macならg77で解く
だった気がします。

Fortranは科学技術計算に特化していますが、文法は分っても、「自分が書きたいプログラム」の書き方の参考になるものは少なかった気がします。77だとそんなこと気にすることもないですが・・・。common文は使わない、とかimplicit realは使わないくらいでしょうか。

当時、90/95の書籍はあったのですが、流体の計算をする場合の書き方みたいなものを提示してある書籍かつインストール方法やgnuplotを使った可視化まで書いてある本はあまりなかった記憶があります。この本は自習用として良いのではないかとすぐ購入した覚えがあります。

結果、モジュールやサブルーチンの使い方、プログラム規模に応じた使い分け、浅水流方程式の特性曲線法を使ったプログラムまで掲載されており、手垢塗れになりました。付箋もいっぱい貼ってあります。

77はGOTO文を使わないとLOOPを抜け出せない、文番号をつかったFORMAT文、1行72文字まで等の制限がありますが、90/95は自由形式になり、それなりに使いやすくなってます。配列も使いやすくなってます。

Fortranを使う人間は、あまりプログラムがメインでない技術者が多く、文法とか書き方にこだわりすぎても仕方ありません。Cだとメモリに気を使わないとダメですし。
その点、Fortranはちょうど良いかと思います。フリーのライブラリも多くありますし。

フリーのコンパイラならgfortran、商用ならインテルコンパイラやPGIコンパイラがあります。Linuxかつ非商用目的なら、インテルコンパイラはフリーで使えます。

あとFortran一番の特徴は計算が速いことです。これが利用される目的として一番大きいですかね。時間の数値積分のための繰り返し回数が多くなると、計算終わるまでに時間がかかります。計算が速いのは重要です。

関数プログラミング

2014-02-10T10:18:00.002+09:00

プログラミングの書き方はいくつかあると思いますが、
私のこれまでの感覚だと、手続き型とオブジェクト指向型がメインだったように思います。

最近の書籍、雑誌、ネットでは「関数プログラミング」というキーワードが多くあるとおもいます。

私が関数プログラミングに触れ始めたのはEmacsからです。Emacsを使いこなそうとすると、Elispを使わなければなりません。いや、使わなくてもEmacsは使えるんですが、これがないなら他のエディタを使った方が良いです。

Elispを見て、それまでの命令型と全く違って面食らいました・・・。
かっこまみれは良いとして、関数型の思想を勉強してないので、頭が大混乱・・・。

そこから関数型言語としてCommonLispとか、Schemeとかあるんだなと。
ハッカーと画家も買って読みました。読み物として面白いです。
が、その威力が未だに見えてきません。変数書き換え不可！なぜ？

ある程度大規模な科学技術計算（私は流体解析メイン）はFortran、C/C++で書かれることが多く、空間離散点が数百万、時間積分も長時間やる場合は、Intel、PGIのコンパイラが多いようです。

中規模であれば、javaやC#なんかも使われているようです。

小規模であったり、統計解析や定常計算等であれば、Excel、Python、Ruby、R等が使われることが多いようです。

さて、自分のフィールドである科学技術計算に関数プログラミングを使えないかと考えたのですが、なかなか難しいです・・・。

その理由の第一は関数プログラミングを十分理解出来ていないことです。これは本読んで、コーディングすれば少しずつ解決しそうです。

第二に関数プログラミングが適用されている他のフィールドで求められている機能について理解出来ていないことです。これはなかなか難しそうですが、第一の理由を解決する過程で何か答えの断片を得られそうです。

その上で、第一、第二を踏まえて科学技術計算にどのように適用するのかを考えようと思います。

一つの機能としては、並列化への適用は可能性がありそうです。科学技術計算をやる上で、「解析時間の短縮」は大きなウェートです。
可読性向上については、厳しいような気がします。計算屋は命令型しか勉強しないですし、オブジェクト指向すら知らない人が多いです。実際、オブジェクト指向で書くと、遅くなることが多いです。
あとは、保守管理、過去の遺産の利用等々、いくつかの側面についても検討する必要があると思います。

ちなみに、関数プログラミングについては、Twitterで、おそらくその道のプロ（NY先生）であろう方からの返信があって、考え始めました。資料もご提示いただき、少し頭が晴れてきました。このページを見ているとは思いませんが、お礼申し上げます。ありがとうございます。

また、考えがまとまったら書こうと思います。

KKスキーム

2014-02-09T01:39:00.003+09:00

河村・桑原スキームについてです。

1984年のLESに関する論文だったかと思います。

数値解析のうち、移流項の解法は以下の流れで勉強しました（私の場合）。

・中央差分（必ず発散）
・1次精度風上差分（解けるけど解が鈍る）
→中央差分の方が精度が良さそうだが解けない。なぜなら移流項（対流項）は移流速度に乗って情報が上流から下流に伝搬するからである。1次精度風上差分の場合、意図しない所で、拡散的な項（数値粘性）が入っているため解が鈍る。

・LaxWendroff法（2次精度だが振動が入る）
・MacCormack法（2次精度だが振動が入る）
→1次精度風上差分では数値粘性が入っていたため、安定的に解けた。同様の考え方で人工粘性を入れることで、振動を小さくする。

・KKスキーム
→3次精度を実現するために4次精度の中央差分＋4次の拡散項を入れる。振動は残るが、河村先生の論文以降、実現象解明へ適用されてきた実績がある。

以降、QUICK、TVD（MUSCL系、Non-MUSCL系）、CIP（オリジナル、有理関数補間、CSL系）の勉強をしました。個人的にはCIP-CSLR1が好きですが、既存のコードの書き換えは結構大変です。マルチモーメントなので・・・。
その後は、小松先生の6Point Scheme、牛島先生のQSIスキーム、ENO、WENO等もかじってます。

今回、気が向いて、KKスキームについてちょっと勉強し直しました。
ベンチマークで使われる2D-Cavity Flowに適用しようと考えてます。
非圧縮性流体の解析なので、MAC系の解法を使って、移流項にKKスキームをと・・。

KKスキームは非保存系のスキームでした。ということで、ある論文の方法を使って保存系に変形しました。線形移流方程式に対して、保存系、非保存系で解いた結果、同様の結果を得られたので、適用は出来そうです。

うまくいったら、結果を載せようと思います。

xojo

2014-02-04T21:30:00.001+09:00

xojoの2月の課題があったのでやってみた。
全然きれいじゃないし分ってないことを認識しているが、
とりあえずやってみみた。

下記ページを参考にさせていただいた。

http://www.geocities.co.jp/SiliconValley/2413/tips.html#16

なんかプリントスクリーンが巧くいかない・・・。

河川のシミュレーション!

2014-01-24T08:17:00.001+09:00

河村先生の書籍です。
現在は絶版になっております。

大学の時に購入しました。
解析の基本的な説明の後、1次元、2次元、3次元の解析について記述されています。

偏微分方程式については移流拡散方程式、ポアソン方程式の解法を、
また一般座標系の説明があります。

河川に焦点を当てた（浅水流方程式に焦点を当てた）書籍は他にないような気がします。

サンプルプログラムもありますが、実用性は低いように思います。
また、流砂量公式については記述がありません。

ただし、勉強し始めた人にとっては有用な本の一つであることは間違いないと感じました。最近は大学レベルでも1からプログラムを書くことは少なくなっているようですし、
あくまでソフロウェアのユーザーであっても、知っておいた方が良いように思います。

個人的にはKKスキームについて記述されていないのが、ちょっと残念です。
日本人が作った数少ない移流項の解法ですからね。
（河村先生のその他の本には記述されています。河川だから省いたのかなと邪推しています）

日本の技術書は薄い本が多いように感じます。
ゼロから勉強しようにも、技術のバックグラウンドや基礎方程式について記述されているだけの本が多いです。基本的には基礎方程式は複雑なことが多く、複雑な方程式を解くには数値計算に頼らざるを得ないことが多いです。そのためには数値解析の基本的事項およびそれを実用的に使う方法まで書くには厚くすべきだと思ってます。
その点、海外の書籍は厚く、読むのが大変ですが、独学に向いています。
日本でもそのような流れにしていただきたいです。

ちなみに河村先生は多くの書籍を出されています。数値計算の基本的な所を記述した教科書になるものが多いです。私はファンなので買ってしまいます・・・。
エクセルVBAを使った解析について書かれているものもありますので、興味のある人は是非！

履歴

2014-01-24T07:00:00.001+09:00

年明け一発目です。
せっかくなので、これまでのプログラミング履歴をまとめておこうと思います。
順番は適当です。

C/C++

Cは大学のときに最初に学んだ言語。その後、数値解析に目覚めて、CIP法を書いてみて、数式を展開して、思いどおりの結果が出てハマる。研究室に所属して、修論ではC++となるが、当時はあまりオブジェクト指向を理解していなかった。社会人になってからはチマチマ勉強して多少理解が進んだ。業務でも扱うこと(VC++)が多く、とりあえず、一番利用していると思う。

Fortran

科学技術計算をやるならFortranは外せない。もちろんゼロから書くならC/C++でも良いが、如何せん、研究室の遺産を利用しようとすると、書けなくても読めなければならなかった。まあ、C/C++に比べれば大したことはない。研究室では77しかなかったが、社会人になってから京都大学の牛島先生の本が出て、90/95を勉強した。当時は配列の動的確保やintent文が使えるなという印象だった。牛島先生の本は90/95の文法だけでなく、科学技術計算のプログラムの書き方について記述してあるのが優れていると思う。

python

初めて存在を知ったのは、研究室に入ってから。共有のパソコンにヘビのアイコンがあったのを記憶している。しばらく忘れていたが、社会人になってみんpyをちょくちょく読みながら少しずつ勉強した。高速な計算でなければ使いやすいと思う。数学者が作成した言語というのも好感度アップ。matplotlibによるグラフ作成や、ベクトルやコンター図も簡単にできるのが良い。また、マルチパラダイムなのも良い。特にFortranユーザーにも使いやすいと思う。海外研修時に簡単な温度応力解析のプログラムをこいつで書いた。

ruby

日本人が作ったとのことで、たのしいRubyを手に取ったのが始まり。pythonに比べて科学技術計算用途のmoduleは少ないと思う。ただし、オブジェクト指向の勉強としては優れていると思う。Railsは使ったことがない・・・。ごめんなさい。データ整理なら特に問題なく使える。業務で使ったことはない。ただし、日本人ならruby使って発展させた方が良い気がする。黒魔術も使えるし。

perl

仕事で海外研修時にO'Reillyの例のやつ（初めてのPerl）を読んでいた。出張中はWindowsしか触れなかったのでPadreを入れていたと記憶している。言語学者が作っているのであまりピンとこなかった。ただし、上記書籍は読んだ方が良いと思う。普通に面白い。業務では科学技術計算用Linuxの設定のついでにLAMPを入れた時に多少触ったくらい。python、ruby、perlについては用途が被るが好みで選択すれば良いと思う。どれが優れているとかは感じなかった。

xojo

real basicだったのxojoになった。勉強会も開かれているようだ。もちろん会社では一切使わない。移流方程式を解いて、図化する程度のプログラムは書いてみた。これも海外出張時に多少勉強した。日本語の情報は少ないが、英語の情報はたくさんあるし、リファレンスやサンプルはしっかりしている。ただ、サードパーティなので一抹の不安はある。Linux、Mac、Windowsの環境で使えるし、ネイティブコードを吐き出せるのは良い所。

いや、プログラム言語じゃない！と言われると思いますが、ごめんなさい。会社ではデータ整理やちょっとした統計解析なんかはVBAでやることが多いです。VBAが優れている点は、Excelがインストールしてあれば細かい設定は関係無しにみんな使える所です。通信関係は触ったことはないです。

VB.NET

メールを読み込んで添付ファイルを抽出するプログラムを業務で書いた。好んで使った訳でなく、遺産がVBだったから使ったまで。文法については「Must Inherits」とかカッコいいと思った。6.0から.NETで強烈に変化してしまったため、ユーザーは減ってるらしい。私の周り？そもそもプログラム書くことに情熱を持ってるやつなんていないっす。

java

javaによるCIP法の本が出てちょっとだけ勉強した。解析結果を表示しながら解析を進めるというもので、結構面白かった。ただ、Fortran的な書き方しかしてないので、javaを勉強した気はしなかった。ほんとはちゃんと勉強したい。Androidアプリも弄ることが出来そうだし。。。

ELisp

Emacsを使い始めると必然的に使わなきゃいけない。今だに使えるようになったとは思わないが、Emacsは楽しい。C/C++から入った私としては、文法に面食らった。何かに使うわけでなくEmacsのカスタマイズのみである。なおWindows環境でのEmacsライクなエディタにxyzzyがあるがこちらはCommon Lispでカスタマイズする。

読んだ本はプログラミング言語AWK。スクリプト言語はsh、sed、awk等が基本的な所にあると思う。利用としてはテキスト処理しかないが、結構簡単に書けるしパイプラインで繋いで処理していく過程は楽しい。完全な主観だが、ワンラインならはやりのスクリプトよりこっちの方が良いかもしれない。

C++

2013-12-18T00:38:00.001+09:00

C++はマルチパラダイム言語ですが、
最近、再勉強をしています。

楽しいです。とりあえず、微分方程式を解くためのアルゴリズムを中心に書いていますが、本当に色々な書き方があるんだなと。

本も数冊購入しました。一冊は本当に初学者の本です。なぜなら、C++の全体のイメージを思い出すことが必要だと思ったからです。新しい発見でしたが、比較的初学者の本を読んで、ちょっと説明が違うな、と認識できるところがありました。少しは成長していたんだなと。

別の一冊はSTL中心です。STLは本当にコードが短くなるだけでなく、デバッグすべき項目を減らせるなという印象です。世の中の、どの分野で、どの程度使われているかわかりませんが・・・。加えてBoostの知識も入れておこうと思います。

最後の一冊はWindowsに特化したMFC中心の本です。
今さら感が強いですが・・・。
VC++はC++とは全然違うものというのが再認識されました・・・。

ライブラリを含めると、C++という言語はあまりにも広すぎます。
プラットフォームに依存したVC++だと別物です。

なかなか学習の仕方が難しいですが、今のところ、やりたいことは微分方程式を解くことなので、いろいろな書き方を試したいと思います。

全然書けてない・・・。

2013-12-12T17:51:00.002+09:00

ようやく落ち着いたので、書き始めようと思います。

次の会社に向けて、C++の復習と統計の勉強を少しずつしています。

・C++について
プログラミング言語の学習は文法をマスターすることではありません。
その言語がどのような思想で設計され、与えられる課題に対して、どのようにアプローチすべきか？というような「思考方法」を学ぶことだと捉えています。

C++のようないろんな書き方がある場合、つまりは思考方法が複数あることを意味しています。そのような混沌としたような言語は非常に難しく感じます。
ただし、それが非常に面白いです。

結局、C言語の思考方法に振り返ることになってます（笑）
でもコーディングする人間が忘れてはならないことは、
「このコードが何をするために使われるか？」
だと思います。

・統計
「統計学は最・・・」とかいう本が出てますが、読んでません。
個人的には、統計学で説明できる事象とは何か？に興味があります。
特に時系列分析やスペクトル解析でしょうか。

統計学はあくまで考え方の一つであり、統計学に准ずれば正解、ということではありません。如何せん、現象論は考えていないからです。ただ、現象としては捉えられるけど、メカニズムは不明、という状況において、何らかしら定量的なアプローチをする方法としては有効であるように思います。すべての現象をニュートン力学や微分方程式で表現できる訳ではない（いや、ほとんどは出来るんですが、その場合、時空間解像度を相当確保しなければなりません）ですし、表現できたとしても、実現象を再現できるほど精度がある訳ではありません。

そういった現象と結果を結びつけたり、そのためのツールとして非常に有効な学問だと思います。平均値、相関、標準偏差とかです。

でも、中学や高校で確率の勉強が、工場の品質管理に使われるなんて想像もつきませんでした（笑）。鵜呑みにして勉強した結果、あ、そういうことか、と納得したのが大学です。

さらに言えば、水文学なんかもそうです。母数推定や関数へのフィッティング等。。。

怖いのは、ツールを使えばなんか答えがでる、ということです。
勉強していない人間がテキトーなことをして、定量的に評価することほど怖いことはありません。

統計学にはそういう恐ろしさも孕んでます。

RCIP

2013-10-20T19:07:00.000+09:00

CIPは非常に優秀な移流方程式に対する解法ですが、不連続部に若干のオーバーシュート、アンダーシュートが残ります。

これは補間関数に3次関数を選択したためであり、これを避けるために、別の関数を使えば良いじゃん！というのが、RCIP(Rational CIP)基本的な考え方です。有理関数を使えばオーバーシュートアンダーシュートを抑えられます。

ソースを晒しておきます。まあ、汚いです。

#advection #################################################

# RCIP ########### ########################################

############################################################

#import#####################################################

from pylab import *

#set parameters ############################################

in_ff = open('test.txt')

imx=100

dx=1.0

dt=0.2

u=0.2

alpha = 1.0

eps = 1.0e-15

print "-"*40

print "u =",u

print "dx =",dx

print "dt =",dt

print "-"*40

############################################################

#initial condition ################################

f= [0]

fn= [0]

df = [0]

dfn = [0]

fini=[0]

fl= [0]

fr= [0]

x= [0]

abm= [0]

xl= [0]

f_tilde = [0]

fex = [0]

i=1

for x in in_ff:

print i

a = x.split(" ")

for k in range(0,1):

print a[k]

xl=xl+[float(a[0])]

f=f+[float(a[1])]

fini=fini+[float(a[1])]

fn=fn+[float(a[1])]

fl=fl+[float(a[1])]

fr=fr+[float(a[1])]

f_tilde = f_tilde + [0]

fex = fex + [0]

df = df + [0]

dfn = dfn + [0]

i=i+1

imx=i

adm=[0.0]

for i in range(1,imx):

adm=adm+[0.0]

print "*"*40

print "xl,f"

for i in range(0,imx-1):

print float(xl[i])," ",float(f[i])

# set parameter ####################################

for n in range(0,1000):

###############################

for i in range(1,imx-1):

if u>=0.0:

ip = i - 1

dd = -dx

else:

ip = i + 1

dd = dx

xi = -u*dt

S = (f[ip]-f[i])/dd

if S - f[ip] >0.0:

isgn = 1.0

else:

isgn = -1.0

TS = f[ip]-S

BB = ((abs(S-df[i])+eps)/(abs(S-df[ip])+eps)-1.0)/dd

a3 = (df[i]-S+(df[ip]-S)*(1.0+alpha*BB*dd))/(dd*dd)

a2 = S*alpha*BB+(S-df[i])/dd-a3*dd

a1 = df[i] +f[i]*alpha * BB

a0 = f[i]

fn[i] = (a3*xi*xi*xi+a2*xi*xi+a1*xi+a0)/(1.0+alpha*BB*xi)

dfn[i] = (3.0*a3*xi*xi+2.0*a2*xi+a1)/(1.0+alpha*BB*xi)-(alpha*BB)*(fn[i])/(1.0+alpha*BB*xi)

#shift#########################################

for i in range(1,imx-1):

f[i]=fn[i]

df[i]=dfn[i]

print fn

###################################################

# Exact ###########################################

for i in range(0,imx-1):

if i>=50 and i<=60:

fex[i] = 1.0

else:

fex[i] = 0.0

#output############################################

ff=open("output.dat","w")

for i in range(0,imx-1):

ff.write(str(xl[i]))

ff.write(" ")

ff.write(str(f[i]))

ff.write("\n")

ff.close()

#output

#for i in range(0,imx-1):

# print f[i]

print "-"*100

print xl

print f

#####################graph#########################

y1=f

y2=fini

y3=fex

plot(xl,y1,label = 'RCIP')

plot(xl,y2,label = 'INITIAL')

plot(xl,y3,label = 'EXACT')

xlabel('x')

ylabel('f')

title('RCIP')

axis([-0,100.0,-0.2,1.2])

legend()

grid(True)

savefig('RCIP')

show()

###################################################

等間隔WENO

2013-10-20T13:07:00.003+09:00

WENOは格子間をラグランジュ補間して、スムーズインジケーターで修正する手法らしいです。ただ、あくまで私の理解であり、正しい認識を提供するものではありません。

河川の流れを対象とし、RANSを用いた計算だと、WENOのような高級な手法を用いることはありませんが、勉強がてらやってみました。

汚いですが、ソースも晒しておきます。たぶん間違いがあると思いますので、あくまで参考ということで。

あとは不等間隔の場合、係数が異なるように思いますので、座標変換して解くような手法（格子間隔を一定として変換）でない限り、修正が必要かと思います。QUICKと同様ですね。

#advection #################################################
# WENO ########################################
############################################################
#import#####################################################
from pylab import *
def calflux1(u,f,flm,flp):
alpha = 0.0
for i in range(1,imx-1):
alpha_tmp = u
if abs(alpha_tmp) > abs(alpha):
alpha = alpha_tmp
for i in range(1,imx-1):
flm[i] = 0.5*(f[i]*u - alpha*f[i])
flp[i] = 0.5*(f[i]*u + alpha*f[i])
flm[0] = 0.0
flm[imx-1] = 0.0
flp[0] = 0.0
flp[imx-1] = 0.0
def calflux2(i,omegam,omegap,flm,fip,fl2m,fl2p):
print i
fltmpp = [0]
fltmpm = [0]
for k in range(0,3):
fltmpp = fltmpp + [0]
fltmpm = fltmpm + [0]

fltmpp[0] = 1./3.*flp[i-2] - 7./6.*flp[i-1] + 11./6.*flp[i ]
fltmpp[1] = -1./6.*flp[i-1] + 5./6.*flp[i ] + 1./3.*flp[i+1]
fltmpp[2] = 1./3.*flp[i ] + 5./6.*flp[i+1] - 1./6.*flp[i+2]
fl2p[i] = 0.0
for k in range(0,3):
fl2p[i] = fl2p[i] + omegap[k]*fltmpp[k]

fltmpm[1] = -1./6. * flm[i-1] + 5./6. * flm[i ] + 2./6. * flm[i+1]
fltmpm[2] = 2./6. * flm[i ] + 5./6. * flm[i+1] - 1./6. * flm[i+2]
fltmpm[3] = 11./6. * flm[i+1] - 7./6. * flm[i+2] + 2./6. * flm[i+3]
fl2m[i] = 0.0
for k in range(0,3):
fl2m[i] = fl2m[i] + omegam[k]*fltmpm[k]
in_ff = open("test.txt")
dx=1.0
dt=0.2
u=0.2
c1 = 1./10.
c2 = 6./10.
c3 = 3./10.
eps = 1.0e-10
print "-"*40
print "u =",u
print "dx =",dx
print "dt =",dt
print "-"*40
############################################################
#initial condition ################################
f= [0]
fn= [0]
fini=[0]
fl= [0]
fr= [0]
fex = [0]
x= [0]
abm= [0]
xl= [0]
flm = [0]
flp = [0]
flux = [0]
flux1 =[0]
flux2 =[0]
flux3 =[0]
fl2m = [0]
fl2p = [0]
i=1
for x in in_ff:
print i
a = x.split(" ")
for k in range(0,1):
print a[k]
xl=xl+[float(a[0])]
f=f+[float(a[1])]
fini=fini+[float(a[1])]
fn=fn+[float(a[1])]
fl=fl+[float(a[1])]
fr=fr+[float(a[1])]
i=i+1
fex = fex+[float(a[1])]
flux = flux + [0]
flux1 = flux1 + [0]
flux2 = flux2 + [0]
flux3 = flux3 + [0]
flm = flm + [0]
flp = flp + [0]
fl2m = fl2m + [0]
fl2p = fl2p + [0]

imx=i
adm=[0.0]
alpham = [0]
alphap = [0]
omegam = [0]
omegap = [0]

for k in range(0,3):
alpham = alpham + [0]
alphap = alphap + [0]
omegam = omegam + [0]
omegap = omegap + [0]
print "*"*40
print "xl,f"
for i in range(0,imx-1):
print float(xl[i])," ",float(f[i])
# set parameter ####################################
# 1st order upwind ####################################
for n in range(0,1000):
calflux1(u,f,flm,flp)

# WENO
for i in range(3,imx-4):

isp1 = 13./12.*(flp[i-2]-2.0*flp[i-1]+flp[i])*(flp[i-2]-2.0*flp[i-1]+flp[i])+1.0/4.0*(flp[i-2]-4.0*flp[i-1]+3.0*flp[i])*(flp[i-2]-4.0*flp[i-1]+3.0*flp[i])
isp2 = 13./12.*(flp[i-1]-2.0*flp[i]+flp[i+1])*(flp[i-1]-2.0*flp[i]+flp[i+1])+1.0/4.0*(flp[i-1]-flp[i+1])*(flp[i-1]-flp[i+1])
isp3 = 13./12.*(flp[i]-2.*flp[i+1]+flp[i+2])*(flp[i]-2.*flp[i+1]+flp[i+2])+1./4.0*(3.*flp[i]-4.*flp[i+1]+flp[i+2])*(3.*flp[i]-4.*flp[i+1]+flp[i+2])

ism1 = 13./12.*(flm[i+1]-2.0*flm[i+2]+flm[i+3])*(flm[i+1]-2.0*flm[i+2]+flm[i+3])+1.0/4.0*(3.0*flm[i+1]-4.0*flm[i+2]+flm[i+3])*(3.0*flm[i+1]-4.0*flm[i+2]+flm[i+3])
ism2 = 13./12.*(flm[i]-2.0*flm[i+1]+flm[i+2])*(flm[i]-2.0*flm[i+1]+flm[i+2])+1.0/4.0*(flm[i]-flm[i+2])*(flm[i]-flm[i+2])
ism3 = 13./12.*(flm[i-1]-2.0*flm[i ]+flm[i+1])*(flm[i-1]-2.0*flm[i]+flm[i+1])+1./4.0*(flm[i-1]-4.*flm[i]+3.*flm[i+1])*(flm[i-1]-4.*flm[i]+3.*flm[i+1])
#
alpham[0] = c1/(eps + ism1)/(eps+ism1)
alpham[1] = c2/(eps + ism2)/(eps+ism2)
alpham[2] = c3/(eps + ism3)/(eps+ism3)
#
alphap[0] = c1/(eps + isp1)/(eps+isp1)
alphap[1] = c2/(eps + isp2)/(eps+isp2)
alphap[2] = c3/(eps + isp3)/(eps+isp3)
#
sigmaalpham = alpham[0] +alpham[1] +alpham[2]
sigmaalphap = alphap[0] +alphap[1] +alphap[2]
#
omegam[0] = alpham[0] / sigmaalpham
omegam[1] = alpham[1] / sigmaalpham
omegam[2] = alpham[2] / sigmaalpham

omegap[0] = alphap[0] / sigmaalphap
omegap[1] = alphap[1] / sigmaalphap
omegap[2] = alphap[2] / sigmaalphap

calflux2(i,omegam,omegap,flm,flp,fl2m,fl2p)
fl2m[2] = 0.0
fl2m[1] = 0.0
fl2m[0] = 0.0
fl2p[2] = 0.0
fl2p[1] = 0.0
fl2p[0] = 0.0

fl2m[imx-1] = 0.0
fl2m[imx-2] = 0.0
fl2m[imx-3] = 0.0
fl2p[imx-1] = 0.0
fl2p[imx-2] = 0.0
fl2p[imx-3] = 0.0
# EULER ##############################
for i in range(3,imx-4):#0~imx-1
fn[i]=f[i]-dt*(fl2p[i]-fl2p[i-1]+fl2m[i]-fl2m[i-1])/(dx)
#shift#########################################
for i in range(2,imx-2):
f[i]=fn[i]
print fn
###################################################
# Exact ###########################################
for i in range(0,imx-1):
if i>=50 and i<=60:
fex[i] = 1.0
else:
fex[i] = 0.0
#output############################################
ff=open("output.dat","w")
for i in range(0,imx-1):
ff.write(str(xl[i]))
ff.write(" ")
ff.write(str(f[i]))
ff.write("\n")
ff.close()
#output
#for i in range(0,imx-1):
# print f[i]
print "-"*100
print xl
print f
#####################graph#########################
y1=f
y2=fini
y3 = fex
plot(xl,y1,label = 'WENO')
plot(xl,y2,label = 'INITIAL')
plot(xl,y3,label = 'EXACT')
xlabel('x')
ylabel('f')
title('WENO')
axis([-0,100.0,-0.2,1.2])
legend()
grid(True)
savefig('WENO')
show()
###################################################

HartenYee

2013-10-16T20:52:00.001+09:00

汚いですし、直したいところはたくさんありますが、
晒すことにします。

HartenYeeの2nd Order Upwindです。
おそらく数値計算ハンドブックにも載ってるような記憶がありますが、
基本的には藤井先生の本をトレースしたものです。
「test.txt」は０〜１００の格子番号iに対して10<i<20でf＝1.0、その他で０としています。

あとは[pylab]を入れていればそのまま描画されるはずです。
ない場合は、import以下を消して最後のグラフ描画部分も消して下さい。

数年前のプログラムに、修正を施しただけですので、
もっと、きれいに書けます（もっとOOPな感じで）が、まあ、残しておきます（笑）

OOPならC++、Rubyで書いてますが、Rubyが一番オブジェクト指向にしっくりきます。

ただ、科学技術計算ならPythonが一番使いやすいと思います。描画は特に。

ちなみにHartenYeeのスキームはあまり利用されていないようですが、
Non-MUSCL系だと好きだったりします。

「中央差分に修正流束を加算する」というスキームの形になってるのが理解しやすいです。

#advection #################################################
# Harten Yee Upwind ########################################
############################################################
#import#####################################################
from pylab import *

# Harten Yee ##############################################
def cal_fai(z):
eps = 0.0001
if abs(z) > eps:
return abs(z)
else:
return (z*z+eps*eps)/(2*eps)
def cal_sigma(z,dt,dx):
return 0.5*(cal_fai(z)-dt/dx*z*z)
def minmod(x,y):
if abs(x)<=abs(y) and x*y >=0.0:
return x
elif abs(x)>abs(y) and x*y >= 0.0:
return y
elif x*y<0:
return 0.0
else:
print "Error minmode2"
#set parameters ############################################
in_ff = open('test.txt')
imx=100
dx = 1.0
dt = 0.2
u = 0.2
print "-"*40
print "u =",u
print "dx =",dx
print "dt =",dt
print "-"*40
############################################################
#initial condition ################################
f= [0]
fn= [0]
fini=[0]
fl= [0]
fr= [0]
x= [0]
abm= [0]
xl= [0]
f_tilde = [0.0]
limiter = [0.0]
fex=[0]
df = [0]
gi = [0]
gamma = [0]
fai = [0]
i=1
for x in in_ff:
print i
a = x.split(" ")
for k in range(0,1):
print a[k]
xl=xl+[float(a[0])]
f=f+[float(a[1])]
fini=fini+[float(a[1])]
fn=fn+[float(a[1])]
fl=fl+[float(a[1])]
fr=fr+[float(a[1])]
f_tilde = f_tilde+[0.0]
fex = fex+[float(a[1])]
limiter = limiter+[0.0]
df = df + [0.0]
gi = gi +[0.0]
gamma = gamma + [0.0]
fai = fai + [0.0]
i=i+1
imx=i
print "imx=",imx
print "*"*40
print "xl,f"
for i in range(0,imx):
print i," ",float(xl[i])," ",float(f[i])
# Time integration ################################
for n in range(0,1000):
#LAX-WEDNROFF##############################
for i in range(0,imx-1):
df[i] = f[i+1] - f[i]
for i in range(2,imx-1):
gi[i] = minmod(df[i],df[i-1])
for i in range(3,imx-2):
if df[i]!=0.0:
#print i,df[i],gi[i],gi[i+1]
gamma[i] = cal_sigma(u,dt,dx)*(gi[i+1]-gi[i])/df[i]
else:
gamma[i] = 0.0
for i in range(2,imx-1):
fai[i] = cal_sigma(u,dt,dx)*(gi[i]+gi[i+1])-cal_fai(u+gamma[i])*df[i]
for i in range(1,imx-3): # Numerical Flux
CS = 0.5*(u* f[i+1] + u*f[i] ) # CS:Centered Space
CF = 0.5*fai[i] # CF:Correct Flux
f_tilde[i] = CS + CF # Numerical FLUX
#print i,limit
for i in range(2,imx-2):#0~imx-1
fn[i] = f[i] - dt/dx*(f_tilde[i] - f_tilde[i-1])
#shift#########################################
for i in range(2,imx-2):
f[i]=fn[i]
print fn
###################################################
# Exact ###########################################
for i in range(0,imx-1):
if i>=50 and i<=60:
fex[i] = 1.0
else:
fex[i] = 0.0
#output############################################
ff=open("output.dat","w")
for i in range(0,imx-1):
ff.write(str(xl[i]))
ff.write(" ")
ff.write(str(f[i]))
ff.write("\n")
ff.close()
#output
#for i in range(0,imx-1):
# print f[i]
print "-"*100
print xl
print f
#####################graph#########################
y1=f
y2=fini
y3=fex
plot(xl,y1,label = 'Harten-Yee')
plot(xl,y2,label = 'INITIAL')
plot(xl,y3,label = 'EXACT')
xlabel('x')
ylabel('f')
title('Hatren-Yee')
axis([-0,100.0,-0.2,1.2])
legend()
grid(True)
savefig('Harten-Yee')
show()
###################################################

MUSCL

2013-10-10T22:51:00.001+09:00

MUSCLをVBAで書いてみました。
Minmod関数でリミッターを掛けてます。

緑色はセル境界の内挿した値です。

河川ならこの程度の精度があれば十分だと解釈してます。
（個人的な意見です）

ただし、移流項に対してのみTVD条件を満たす（システム系で解かずに、運動方程式の移流項に対してのみTVD条件を満たす）場合、他の項による不連続部に対して、リミッターがかかってしまうのではないかと気になっています。

そうすると、熱移動と流れの数値解析や数値流体工学に議論されているハイブリット法が良いのか？と思いますが、アレは拡散項ありきなんですよね。河川だと拡散項が省略されることが多く、生成項が効いてきます。ハイブリッド法だと1次風上のみになってしまいます。

ここら辺の議論については非常に興味深いです。誰も議論しないのだろうか？いや、精度的には意味がないのか？

おそらく河川の流れの空間刻み幅は小さくしてもせいぜい数mなので、高次精度のスキームは必ず必要だと思ってます。

勉強します。

空間２次精度

2013-09-29T22:29:00.002+09:00

空間２次精度についての話です。

河川の流れであれば、
1次精度風上差分では不十分でも、せいぜい２次精度があれば十分なようです。
ダムブレークのような問題を解く必要はあまりないですから。

２次精度だとベースをオリジナルのLaxWendroff法（2STEPではない）を用い、FluxLimiterを施すことで、大局２次精度を維持できます。

CIP法を用いた非保存系のスキームも使われているようですが、既存のプログラムがある場合、書き換えるのが煩雑です。
また、CIP法だとわずかに数値振動が残りますがTVDだと1次精度となって振動を抑えます。
TVDだとMUSCLも使われているようですが、TVD-LaxWendroffは使われることは少ないようです。
精度の高い乱流モデルを後々組み込むことを考えると、MUSCLという選択もわるくないと思います。

以上の話は、ナビエストークス方程式（or レイノルズ方程式）のうちの移流項に対する解法であり、
運動方程式、連続式をシステム系としてとらえた解法ではありません。なぜなら、生成項となる摩擦項や地形勾配項があるためであり、
C-Propertyを満たすことが難しいためです。
アカデミックの世界でも研究されている方はいるようですが、マニアックなまでに研究されている方はいないようです。
すでに解決された問題だと捉えられているのでしょう。
まあ、分っている人は少ないと思いますが・・・。

流体の数値計算法について（２）

2013-09-29T16:45:00.001+09:00

藤井先生の本についての続編です。

一つ面白い記述に

「最近の多くの実用的な計算結果を見る限りでは、一般的には2次精度を維持することが必要といえる」

というのがあります。

乱流の解析では2次精度では不十分な場合がありますが、多くは2次精度程度である程度事足りるとの指摘です。

河川の流れの解析のように、平面2次元流れの解析であれば、この指摘は当てはまると思われます。

平面2次元流れは、どの程度の流れの解像度を得ることが必要かで乱流モデルを選択する必要がありますが、多くの場合は、ブシネスク近似を使って、ゼロ方程式を用いたモデルが多いようです。

①格子形状
格子形状は境界適合座標系を使い、スカラー量とベクトル量を互い違いに配置するスタッガード格子を採用します。

②解析手法
有限体積法を用いることが多いですが、モデルによっては有限差分法が多いです。有限要素法は滅多にありませんが、アカデミックな世界では実施されていることもあります。
また、移流項に対しては、1次風上差分が非常に多いはずです。藤井先生の本を読まれる方は「は？」と思われるかもしれませんが、実用的な解析では、計算が飛ばないことが非常に重要です。
はっきり言って、解析についてあまり知らなくても、計算結果を確実に得ることが求められます。２次精度以上を確保しようとすると、不安定になる可能性が多くなりますので、1次風上が多いのです。
一方、移流項以外は、２次精度の中央差分が多いです。各項の大きさについての議論は別にして、移流項で3次精度以上を確保しようとすると、他の項も同様に3次精度以上を確保する必要があります。
したがって、移流項に対しては大局的空間２次精度程度を確保することが適切のようです。

さて、藤井先生の本に戻りますが、前述の指摘がある程度的を得ています。
分野が違えど、面白い指摘だなと思った次第です。

なお、ブシネスク近似のうちk−εやLESを用いる場合はこの指摘は当てはまりません。
乱流モデルは、小生の間隔では3次精度以上が必要のようです。

書籍

2013-09-16T21:25:00.002+09:00

「プログラミング言語AWK」
AWKの本です。最近、USP Labから出版されたものです。
もちろん著者はAWKの開発者であり、AWK版のK&Rと考えてもらえれば良いとおもいます。

ただ、出版といっても、絶版しては再販されている書籍であり、
需要がある書籍と言えましょう。

読んだ感想ですが、古くさいです。
ただ、K&Rよりも実用的であると感じます。
スクリプト言語であり、ワンライナーによるフィルタをかけるだけであれば、序盤を読めばかなり書けるようになると思います。

まあ、敢えてAWKを使う人はあまりいないと思います。今であれば、Perl、Python、Rubyを使えば比較的簡単にできます。
それがテキスト処理だとしてもです。

ただし、プログラミングの歴史を学ぶ点においては、知っておいて損はないかと。
加えて、「プログラミングが、なぜ、今の形になったのか？」ということを理解する一つの助けとなる気がします。

やはり、周辺知識は独立して使うことが出来るものではないですが、メインの知識と組み合わさって威力を発揮するものだと思います。

等間隔QUICKEST

2013-09-16T02:40:00.002+09:00

# 久々にpython

for i in range(3,imx-3):#0~imx-1
c = u*dt/dx #
if u >0:
fn[i] = f[i] - c/6*(2*f[i+1]+3*f[i]-6*f[i-1]+f[i-2]) \
+ c*c/2*(f[i+1]-2*f[i]+f[i-1]) \
- c*c*c/6 *(f[i+1]-3*f[i]+3*f[i-1]-f[i-2])
else:
fn[i] = f[i] - c/6*(2*f[i-1]+3*f[i]-6*f[i+1]+f[i+2]) \
+ c*c/2*(f[i+1]-2*f[i]+f[i-1]) \
- c*c*c/6 *(f[i-1]-3*f[i]+3*f[i+1]-f[i+2])

等間隔にしか使えないですが・・・。
一応、保存則を使っています。

明日中には不等間隔＋TVD化されたQUICKを晒したいです。

書籍

2013-09-13T06:10:00.001+09:00

藤井先生著「流体力学の数値計算法」について書いておこうと思います。

双曲型の数値計算法（保存性を意識した差分法）がメインです。基礎方程式は圧縮性流体を想定したオイラー方程式です。解析手法はTVD系です。したがって、非圧縮性流体の解析で使われるスタッガード格子やMAC系の解法に関する手法は述べられてません。

ハッキリ言って、これを読んでも実際に解析が出来るようにはなりません。
実用的な解析をするならもっと簡単な本を読んだ方が良いでしょうし、初心者は読まない方が良いかもしれません。

ただし、味わい深い本だと思います。
小生が購入したのは2004年頃ですが、最近ようやく「あ、これのこと言ってたんだ！」と気づかされることが多いです。FDSやFVSの基本的な考え方が書いてある本はこれくらいだと思います。
浅水流方程式への応用はちゃんと鉛筆動かさないとダメですけどね。
正直、まだ、すべてを理解してません（笑）
ただ、「わかんねぇなー」と思いながら、日々少しずつ理解が深まるところに、楽しさを感じてます。

ところで、

「わかりやすいほにゃらら」
「初心者のためのほにゃらら」
「10日でほにゃらら」

等々。分りやすいのでスイスイ読めますが、長くは読めません（否定している訳ではありません）。

結局、簡単に理解できることは、それだけの価値しかありませんし、他人もすぐ理解できます。
苦しんで覚えたことや、時間をかけて学んだことは忘れませんし、他人もすぐには理解できないことです。

良書を見分ける一つの指標として、「古いけど絶版にならない技術書は良い」というのがあります。
ただし、理解に時間がかかるものが多いです。

技術書に限らず、死ぬまでに多くの良書に出会いたいですね。

蛙跳び2

2013-09-04T00:45:00.003+09:00

こんな感じでした。鈍るんです。

LeapFrogは空間積分の右辺をf^{n+1}とf^{n-1}の平均値である0.5*(f^{n+1}とf^{n-1})を使うような良い加減な方法なんであまり好きではないのですが、これ見ると、RungeKuttaを使って計算コスト書けるなら、LeapFrogの方が楽だし、鈍らないし、と判断してしまいますね。何か間違ってることを信じたいです（笑）

蛙跳び

2013-09-03T13:42:00.000+09:00

' Leap Frog Method
dim n as Integer
dim i as Integer
dim f(100) as double
dim fo(100) as double
dim fn(100) as double
dim mojiretsu as String

'Initial Condition
For i = 1 to ixmx
if i > 10 And i<20 then
fo(i) = 1.0
f(i) = 1.0
fn(i) = 1.0
else
fo(i) = 0.0
f(i) = 0.0
fn(i) = 0.0
end if
next i
'Time Marching
for n =1 to itmx
'Leap Frog
for i = 2 to ixmx
fn(i) = fo(i) - 2*u*dt/dx*(f(i) - f(i-1))
next i

'Shift Variables
for i = 2 to ixmx
f(i) = 0.5*(fn(i)+fo(i))
fo(i) = fn(i)
next i
next n

'Output
for i = 1 to ixmx
mojiretsu = Str(i) +"," + Str(f(i))
ListBox2.AddRow(mojiretsu)
next i

4次精度ルンゲクッタと比較すると、2次精度の蛙飛びの方が、鈍りが少ないです。

なぜなんだろう・・・。

ListBox

2013-09-02T22:17:00.000+09:00

//ListBoxのクリア

Listbox1.deleteAllRows

Winで作ったバイナリファイルをMacで初めて開きました。
OS依存の機能を使ってないことはありますが、特に問題なく開きました。

ただし、r2でも日本語入力がおかしい箇所がありました。
まあ、日本語使わなければいいんですけどね（笑）

ホビー用途なら、すごくいい環境ですが、
日本語の解説書が充実しないと広まらない気がします。
勿体ない・・・。

値への変換とか（Xojo）

2013-09-01T06:08:00.000+09:00

FortranやCを使ってると、値の読み込み時点で型指定をするが、
その他の言語だと、すべて文字列で読み込んで、キャストすることが多いみたいです。

Xojoの型変換は下記のとおり。

Str(数値等)：値をstringに変換
Val(文字列)：文字列を値に変換

使わないとすぐ忘れます・・・。

ところで、下記に値のやり取りの記録を書きます。

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1．ListBox1を配置して、表示の場合。

//ListBoxの表示をすべて消す。

ListBox1.deleteAllRows

//ListBox1に文字列Program START!を追加
ListBox1.AddRow("Program START!")

2．TextField1に書き込まれた値を変数に代入

//TextField1に書き込まれた文字列を、別途定義した文字列「strings」に代入
//その後、数値型に変換

Dim strings as string
Dim i as Integer
strings = TextField1.Text
i = Val(strings)
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

早くいろいろ晒したいです。