2013年12月18日水曜日

C++

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

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

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

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

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

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

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

2013年12月12日木曜日

全然書けてない・・・。

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

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

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

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

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

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

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

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

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

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

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

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





2013年10月20日日曜日

RCIP


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

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()
###################################################

2013年10月16日水曜日

HartenYee

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

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

2013年10月10日木曜日

MUSCL


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

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

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

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

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

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

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

勉強します。

2013年9月29日日曜日

空間2次精度


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

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

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

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

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




流体の数値計算法について(2)

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

一つ面白い記述に

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

というのがあります。

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

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

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

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

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

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


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

2013年9月16日月曜日

書籍


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

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

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

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

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

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

等間隔QUICKEST


# 久々に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年9月13日金曜日

書籍

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

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

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

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


ところで、

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

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

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

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

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








2013年9月4日水曜日

蛙跳び2


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


2013年9月3日火曜日

蛙跳び


' 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次精度の蛙飛びの方が、鈍りが少ないです。

なぜなんだろう・・・。

2013年9月2日月曜日

ListBox

//ListBoxのクリア

Listbox1.deleteAllRows

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

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

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






2013年9月1日日曜日

値への変換とか(Xojo)

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)
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

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

2013年8月22日木曜日

C++

プログラムを書く時の設計思想みたいなものって、皆さん考えますか?

マルチパラダイムな言語だとここをしっかりしていないと全くダメな気がしてます。

ただし、C with Classesな書き方だと、書きながら修正とかしていくんですかね?
少し疑問です。書きながら、設計を修正しつつ、また書く。
そういった開発方法もありなのかなと思ってます。
(アジャイルとか勉強してなくて発言してすみません・・・)

virtual classなんかはそのためのものな気がしてます。

プログラムを作成するのは難しく、そして面白いですね。

2013年8月21日水曜日

XOJO

新しくなっていたので、更新したら、無言で落ちなくなりました。

やはり、こまめな更新は大事ですね。

数値解析、テキスト処理、WEBアプリ等々ありますが、
速度は別にしてXOJOは全部いけます。

速さは厳しいですが、可能性大だと思ってます。
GUI周りは扱いやすいです。


2013年8月17日土曜日

WEBアプリ

科学技術計算中心だと、WEBアプリに注目しません。

ただし、単純なデータ処理は、社内でサーバー立ち上げて、サーバーで処理して、クライアントにアウトプットする、なんてことが有効ではないかと思ってます。

基本オフィスは入れてるので、VBAでも良いですが、PerlやPHPをサーバーで動かすと、生産性が向上する気がしてます。

得る知識の使い回しも出来そうですしね。

2013年8月16日金曜日

マルチパラダイム

プログラムの書き方は色々あると思います。

科学技術計算だと「支配方程式を離散化する」ことに特化したFortarnが向いてます。

ただし、その後の維持・管理・保守を考えると、OOPを使え、実行速度が見込めるC++になると思います。

一方、WEBアプリなんてFortranは向いてませんし、そんな機能が備わっていません。

目的だけでなく、規模にもよると思います。
ちょっとした処理をしたい、書きっぱなし、の状況で仰々しく設計することは時間の無駄です。

結局、プログラムの書き方は、状況に応じて使い分けるべきです。
そして、「これだけ覚えてれば大丈夫!」なんてことを思っちゃダメだと考えてます。

継続して、勉強する、創造する、をやっていかないとなぁと思ってます。


2013年8月15日木曜日

滞り

全然、投稿できてないです。

はあ。

Perlの某書籍を購入しました。Kindleの例のアレです。

入門者に取っては取っ付きやすいと思います。

PerlはLinuxやMacではデフォルトで入ってる点、そしてWindowsの環境の整えやすい点が良いと思ってます。伊達に歴史のある言語ではありません。
そして、後続に対する影響も大きいです。C使いなら、違和感はあまりありません。
非常に高級な言語(読み手が考えなければならないという意味で)だと思ってます。

仕事では使いませんが、楽しい言語だと思ってます。

やっぱりプログラムは楽しくなきゃダメですよね。
そして、今まで作られてきた言語はすべからく、何か言語思想があって、それを楽しむ余裕があると良いですね。そのためには苦しいと感じる勉強も必要ですが・・・。


2013年8月9日金曜日

Perl

Perlの用途はCGIが多いのでしょうか?
私はスクリプトは職業柄、テキスト処理程度です。
本題はFortran、C/C++です。速さ第一なので。

テキスト処理、グラフ作成だと、Pythonが一番だったりします。
なので、ほとんどPerlは使いません。Wikiくらいでしょうか。。。

おそらくPerlでもグラフ作成はできるのでしょうが、やろうと思いません。
科学技術計算向きではないです。

一点だけ、コンテキストの解釈、については、非常にレベルが高いなと思ってますし、
Perlのかっこよさかと思います。
まあ、仕事では使いませんけどね。

2013年8月6日火曜日

コロケート格子

ソースを晒せてません。
C++ではあるのですが、Xojoに書き換えてません。

スキームはセミインプリシットを使ってます。
FVSを使えばもうちょっと楽に陰解法を使えると思ってますが・・・。

2013年8月3日土曜日

無題

Linux(Cent OS 6.4)にXojoを入れてみました。

インストールはダウンロードしてダブルクリックするだけという簡単さでした。

LinuxでGUI環境を整えるならxojoが一番良いかと。
Webアプリもいけますし。

WinならVisual Studio、MacならXcodeですが、xojoはすべてのプラットフォームでいけますし、どのプラットフォームでも、他のプラットフォームの実行ファイルを作成できます。

ライセンスがないと実行ファイルの作成は出来ないですが、デバッグ環境でも良く、勉強したいだけなら、フリーのはずです。(間違ってたらすみません)

サードパーティの可能性を感じます。


2013年8月2日金曜日

Ruby

個人的かつ主観的な意見ですが、pythonよりrubyが好きです。
ただし、仕事ではpythonがメインです。
matplotlibがメインですが・・・。

設計思想はrubyが好きですが、
実利を考えるとpythonという選択をしているかっこ悪い状況です。

他力本願は好きではないので、
なんとか出来たら良いんですが、その能力がないです。


2013年8月1日木曜日

superbee

'Xojo
’引数r
'
  dim r1 as double
  dim r2 as double
  dim r3 as double
  dim rmax as double
  r1 = 0.0
  if 2*r>=1 then
    r2 =1.0
  else
    r2 = 2*r
  end if
  if r>2 then
    r3 =2
  else
    r3 = r
  end if
 
  if r1 >r2 then
    rmax = r1
  else
    rmax = r2
  end if
 
  if rmax > r3 then
    return rmax
  else
    return r3
  end if

2013年7月27日土曜日

レガシー

最近、実行速度を気にしているせいか、コンパイル言語を触ることが多いですが、
書いてて、「最近の言語ならもっとスマートに書けるのに」と思います。

でも、FORTRANとかCでエレガントに書くのが、真のプログラマーかなとも思います。

結局プログラムは「どのようなツールで作るか」の選択であり、成果品の本質ではないはずです(ユーザーからしたら、やりたいことがやれれば何でも良い)。

ただプログラマーはプロフェッショナルなので、速さとか、管理とかいろいろ考えた上で選択したいですね。

だから、昔の言語を知っておくのは大事ですし、今の言語を知るのも大事です。
その上で、なぜこのようなやり方を選択したのか、を明確に説明できるようなロジックを持っていたいですね。

2013年7月26日金曜日

無題


Real Basicのファイル出力です。
変数iと配列f、f_exactをストリングで繋いでラインで出力する。
fは計算値、f_exactは厳密解です。

  'ファイル出力用
  fout=GetFolderItem("output.txt")
  txtout=fout.CreateTextFile
  txtout.delimiter =Chr(13)
  txtout.WriteLine "i,f"
  for i =1 to ixmx
    outs = str(i*dx)+","+str(f_exact(i))+","+str(f(i))
    txtout.WriteLine outs
  next i
  txtout.Close
  'ファイル出力用終わり

2013年7月23日火曜日

Vim

vi/vimを使ってます。
しかたなく。
でも覚えると速いですね。プログラムを書く時のインデント周りの設定が出来ればもっと良くなりますね。


2013年7月22日月曜日

C wtith Classes

c++の書き方は色々あると思います。
better C
C with classes
ジェネリック

速度優先だとBetter Cがやりやすいと思います。
最適化が書けやすいみたいです。

でも、保守性も大事です。

このバランスが一番大事だし、解析屋がご飯食べられるのはこの究極の答えが、Case By Caseだからだと思ってます。

今日はDEMでツールを作れそうだったので、構想だけ考えました。
可視化とかがちょっとめんどくさいです。

でも教養として知っておきたいです。

2013年7月21日日曜日

CentOS

インストールが楽ですね。ウィンドウズよりも。

慣れれば、おそらく一番速いですね。
yum install ・・・で全部行けますし。

GUIは使ってないので分りませんが、Qtとか使えるのか試してみたいです。
GUIが必要な訳じゃないですが、作ることができる環境は整えておきたいです。

Apache,PHPはインストールできたのですが、wikiが巧く行きません。
もう少しなんだけどなぁ。

ポメラ2

キーボードはCapsとCtrlを入れ替えられるといいですね。
エディタはEmacsとか使えるといいです。

もう、テキストエディタで突っ走って欲しいです。
乾電池で動くテキストエディタ。
そんな尖った商品、ぜひKING JIMさんには作り続けてほしいです。

2013年7月20日土曜日

ポメラ

ポメラのシャア版を買った。

カレンダー機能は便利だ。コンパイラ入れなくても良いけど、コードのインデントを自動に入れてくれるとプログラム書き放題なのに。

2013年7月18日木曜日

cent os

intelとか入れました。
ゼロから入れるとvi使えないと話にならないですね。

でも非常に楽しい。
オープンスースの方が、可能性が大きいですね。

GPUも入れて、CUDAとかで解析もやりたいです。

CentOSについては備忘録を残したいですね。
ゴールデンウィークにやりっ放しになってるのでなんとかやってみます。

経験者の方から、「これは読め!」みたいなのがあれば教えていただきたいです。

やること

今日はやることがたくさんある。

Linuxの設定。
エミュレータの設定。
物理モデルの確率モデルへの落とし込み。
ネットワークの設定。
シェルによるデータ整理スクリプト作成。
VisualStudioの設定。
業務データの確認。

ちょっときついな。

2013年7月15日月曜日

QUICK

等間隔格子です。
XOJOです。

必要とあればTVD化のソースや不等間隔TVD-QUICKもあります。

'以下ソース

  dim f1 as double
   if u >= 0 then
    f1 = f-u*dt*(3*fdn1+3*f-7*fup1+fup2)/(8*dx)
  else
    f1 = f-u*dt*(3*fup1+3*f-7*fdn1+fdn2)/(8*dx)
  end if
  return f1

2013年7月14日日曜日

Linux

今、CentOSとUbuntuの勉強中です。

RHELとDebianで二つです。

いやー、面白い。ついでにPHPも勉強中。はじめてのPHPとプログラミングPHPの2冊をポイントでゲットしました。ファイル操作関係を巧くできるといいなと思ってます。

情報工学を学んだ方は当然かもしれませんが、
ニュートン力学を基礎とした工学がベースの人間からしたら、コンピュータはツールでしかないんです。ただ、こだわりたいですね。

XOJOは全部同じだから非常に優秀だと思います。
でも、なんかですね。オープンソース使えないとだめだという脅迫観念があります。

皆さんはどうしてるのか気になります。

もっと、解析のソースを公開したいです。

流体が基本ですが、もっと基本の移流、拡散、ラプラス、ポアソン等々、いろんな言語で公開して、みんなで勉強できたら良いなと思ってます。

よろしくお願いいたしますね。


2013年7月7日日曜日

コロケート格子

なんか難しい。。。
巧くイキマセン。。

うまくいったらソース晒すぞ!

無題

なんか、ゾージョーって書くとアクセス増えるので、書かないことにします。

もう限界に近いのか?

①やりたいこと。
②やらなきゃならないこと。
③やることを求められていること。

多分、最低限②、③が社会人の優先順位。
自分は①が一番です。

だから苦しい。

世界を変えるのはゼネラリストじゃない。
一部のスペシャリストだ。
俺はスペシャリストになりたいんだ!

2013年7月4日木曜日

ソース

なるべく公開したいですが、出来てなくてすみません。

勉強したいときに何から手をつけたら良いのか分らないことが多いと思います。

ざっくばらんですが、話題提供を出来れば良いと思ってます。

勉強大好きな人間の成果を世の中に反映し、勉強してる人に対してちゃんとペイできる世の中が本来、日本のあるべき姿な気がします。


dam break

河道地形を踏まえたSource Termの修正がが分ってきた。

後は摩擦ですね。

2013年7月3日水曜日

ダムブレイク4

スタッガードスキームはあまり好きではありません。

なんで、定義場所変えんねん!



2013年7月1日月曜日

ダムブレイク3

いろいろ弄って、運動方程式の記述を間違っていることに気づいた。

修正したけど、急変するところでは、1次風上を使っても、振動が残ります。
ただ、ベースは出来たので、ようやくFDSの人工粘性を入れられます。

Leap Frog、1st_Order_upwind、人工粘性でスーパーロバストにしてやる予定。

2013年6月30日日曜日

ダムブレイク2

膨張波と伸縮波がうまく表現できてないことが分りました。

もう少し、考えてみます。

とりあえずCで書いてますが、いろんな言語で書くにします。

履歴書書かなきゃ・・・。

そろそろ会社に行くかな?

ダムブレイク

DamBreakはオイラーの方程式を基礎方程式として、FDS系の解法が一番良い気がします。コロケート格子で。

スタッガードスキームだと、そのまま離散化すれば良いですが、ショックキャプチャリングがうまくいかなくて、若干の数値振動が残ります。

シビルエンジニアリングの世界ではほとんど1次風上が多いですね。
正直気に食わないですが。。。

FDS系の解法はXojoで書いてみることにします。そのうち、ソース晒しますね。
今のところbetter C的なC++で書いてます。

もう少し考えてみます。


2013年6月27日木曜日

高速化と可読性

高速化はアルゴリズムの修正を行ったり、基礎方程式からは多少外れます。

でも可読性は若干のオーバーヘッドはあるものの、維持管理は楽な気がします。

科学技術計算は特殊な村の気がしてきました。

2013年6月26日水曜日

C

久々にCを書いてみた。
Cの頭に切り替わらない・・・。

Fortranばかり弄ってるから、こんなことになる。
ポインタがあると、自由度がありますね。

数値解析好きとしてはProgram言語なんて2の次のはずですが、
やはりこだわりたいですね。

2013年6月24日月曜日

C言語

0から作るならば、Cを選択する意味はあまりないと思います。しかし、これまでの財産があり、幸か不幸か、FortramとC/C++は実行速度が圧倒的です。

科学技術計算は、むちゃくちゃな格子数、条件で結果を出さなければなりません。なので、Fortran、C/C++が圧倒的です。

ソースコードの保守を考えると手続き型のC/C++の方が有利だと思いますが、C++のクラス作成して、可読性を高めるためのオーバーヘッドとの兼ね合いかと思います。

ただし、一流の解析屋はすべて出来てしかるべきだと思ってます。

解析は構成則とガバニングイクエーションの基本を踏まえてやりたいですね。
モデラーの醍醐味はメッシュ内の挙動をどのように設定するか?
これにつきるようがします。

2013年6月23日日曜日

xojo


 minmode(r)

  if r>=1 then
    r =1.0
  end if
  if r>0 then
    return r
  else
    return 0.0
  end if

vanLeer(r)
return (r+abs(r))/(1+abs(r))

vanLeerは分母が0になっても構わないので使いやすいと思われます。
その他の制限関数は海外のWikiに載ってます。

2013年6月22日土曜日

名称未設定

TVD法のうち、流束制限関数に関して、ソースを晒したいが、まとまってないので、躊躇してます。

そのうちに。QUICK法とかLaxWendroff法とか、今日ではあまり使わないかもしれませんが、非常に勉強になります。

スキームはこだわりを持って使いたいものですね。
「解ければ何でも良いじゃん!」
はちょっと、バカすぎます。
その場合、1次精度のみでロバスト性を求めるだけで良くて、ぼーっと格子点が増やせるハードが出てくるのを待っているだけです。

CIP法(CSL系)が一番好きですが、既存のコードに対して、フラクショナルステップを使って、移流項のみ改善したい場合は、既存のオイラリアンのスキームの方が良いみたいですね。

ゼロから書くなら、CIP法+その他の項があれば完全陰解法にします。
ロバストかつ高精度。


2013年6月21日金曜日

Xojo先ほどの続き

  とりあえず、Actionに下記ソースを加えると、ファイルへ出力できる

 'ファイル出力用
  fout=GetFolderItem("output.txt")
  txtout=fout.CreateTextFile
  txtout.delimiter =Chr(13)
  for i =1 to ixmx
    txtout.WriteLine str(f(i))
  next i
  txtout.Close

LaxWendroff法とXOJO

  Dim s as string
  Dim i as integer
  dim f() as double
  dim fn() as double
  dim n As  Integer
  dim fout as FolderItem
  Dim txtout as TextOutputStream
 
  s = "Program START!"
  lbResult.AddRow(s)
  Redim f(ixmx)
  Redim fn(ixmx)
 
  For i = 1 to ixmx
    f(i) = 0.0
    fn(i) = 0.0
    If i>10 And i<20 Then
      f(i) = 1.0
      fn(i) =1.0
    end if
    s =str(i) +","+ str(f(i))
    lbResult.AddRow(s)
  next i
 
  For n=1 to 100
    for i =2 to Ixmx-1
      '1次風上
      fn(i) = f(i) - u*dt/dx *(f(i) - f(i-1))
      'オリジナルLaxWendroff法
      fn(i) = f(i) -u*dt/(2.0*dx)*(f(i+1)-f(i-1))+0.5*u*u*dt*dt/dx/dx*(f(i+1)-2*f(i)+f(i-1))
    next i
    for i = 2 to ixmx
      f(i) = fn(i)
    next i
  next n
  For i = 1 to ixmx
    s = str(i) + ","+str(f(i))
    lbResult.AddRow(s)
  next i
 
 
 
 

Elisp

自分用にElispをとりまとめています。
Emacsのorgモードを使えば、かなり楽にHTMLで残せます。

ElispはEmacsカスタマイズ用に勉強するつもりでしたが、
なんか今までのプログラム言語に対する概念が変わりつつあります。

テキスト処理、グラフ化、メール処理、WEBの情報収集等々、日常使いでどこまで出来るか手探りですが、楽しんでやって行こうかと。

何よりフリーで使えます。

xyzzyと同様になにか貢献できればいいですね。

2013年6月20日木曜日

Lisp

Emacsのorg-mode。

かなり使えます。
open-junk-fileはElispの勉強になります。

Gaucheの本も買いました。

GUIまで出来れば良いのですが・・・。

言語なんて所詮ツールです。
でも普通のやつの上を行くには、多くのツールに触れた上で、選択したいものです。

ただし、高速の科学技術計算はFortran、C/C++の二択です。
これはWebやってる人には理解できないかもしれませんが・・・。

良い情報を御待ちしております。
小生の勉強の中で御教えできることは、コードを晒したいと思います。

よろしく御願いいたします。

以上

2013年6月19日水曜日

Lisp

ハッカーと画家を読んでます。

面白い。非常に面白い。

もちろん、疑問だったりする所はあるのですが、
面白いんです。

何事も極めると、社会的に価値があると実感します。

あとLispの設計思想みたいなものにも多少触れられると思います。

「母国語はなんですか?」
「Lispです。」

こう答えられたらカッコいいですね。

プログラミング言語は世界共通ですから。

2013年6月18日火曜日

advection

基本的に

1 中央差分に対して修正する(人工粘性のイメージ)。
2 1次風上に対して修正する(1次精度から高次精度にする)。

の二つの考え方かと思います。

どちらでも構わないですが、海外の書籍だと2が多い気がしてます。
一方、日本は1が多いかなという印象です。

高速流の場合、拡散項が相対的に小さいので、移流項のロバスト性は重要かと。
空間に高次精度を用いる場合、時間積分の精度を上げないと不安定になる場合もあります。





2013年6月17日月曜日

Lisp

いろんな方言があるので、どれが一番良いのか。。。

でも優秀なプログラマーはLispを知っているように思います。

たのしくプログラミングしたいです。
だから色々手を出してるような気がします。

ただし、Fortran、C、Lispは圧倒的な気がします。

プログラミングの思想とか、考え方とかはこれらをやれば圧倒的に視野が広がりそうです。

Lispの使い道が今のところEmacsのElispくらいしかない。

科学技術計算は某コンパイラが圧倒的なのでFortran,C/C++の2択しかないですが、
プロトタイプだったり、テキスト処理なんかはLispももっと使える気がします。

あと、Intel Compilerをコマンドプロンプトから楽にコンパイルできるよう、プロンプトを修正しました。短いコードなら非常に早くコンパイルできるとも居ます。
いかせんせん、Visual Studioが重すぎる時はぜひ使いたいと思います。

Schemaを使ってみるか!
xyzzyも使えるので弄ってみることにします。
惜しむらくはWindowsのみとかは勘弁して欲しいですね。

2013年6月16日日曜日

TVD

自由表面流れの乱流モデルを書いてみたが、TVD系のスキームだときれいに乱れが出ない。

やはりCIP、WENO系か(笑)

というのはおいといて、1次精度の風上差分だと駄目だということがよくわかりました。
仕事柄、1次精度でも十分であり、それ以上にロバストであることが重要であることが多いので、2次精度以上でなければならないことが新鮮でした。

まあ、コーディングしたことがない人間が解析やっちゃ駄目ですね。
だいたい「ロバスト」にいろんな意味を入れすぎてます。「猿でもまわせる」のがロバストではないですよ。


2013年6月15日土曜日

QUICK

Quadratic Upstream Interpolation for Kinematics
最近はあまり使われていないようですが、結構好きなスキームです。

TVD化するときはUMIST関数を使いました。
乱流モデル書いてみると、1次風上との違いがよくわかります。

QUICKにはQUICKEST、Ultimate QUICKESTなんかもありますが、
せいぜい2次精度レベルで良いような解析しかしないことが多いです。

2013年6月14日金曜日

elisp

Emacsの設定をもっと弄りたいのでElispを勉強することにする。
Cとはかなり異なるので、面食らうが、面白い。
古くない言語だ。
;;; 変数aに1をセット
(setq a 1) ; => 1

2013年6月12日水曜日

Lisp

Lispで移流方程式を書こうと思い立った。
リストの作成までは簡単そうだが、リストのいじりがうまく行かない。

でも、すごく勉強になります。
ElispでEmacsを弄りたかっただけなのに(笑)



2013年6月11日火曜日

xojo

皆さんxojoしか興味ないのか・・・。

がんばってソースを晒します。

今、Lispなんですけどね。注目は。

MsgBox "Hello,World !"

お決まりですが、必ずチェックしちゃいますね。

Ruby

ご飯食べるのにFortranじゃ無理でしょ。
とりあえずC/C+。
スクリプトにRuby。
Emacsカスタマイズ用にELISP。

行ったり来たりしながら少しずつ前に進めれば良い。

2013年6月10日月曜日

SLIME

何故かインストールがうまく行かん。

そもそも、Linux,Mac,Windowsで同様の環境を整えようとする所が間違ってるかもしれん。

悔しいなあ。じっくりやるか。

2013年6月9日日曜日

解析スキーム

各分野において、解析スキームの種類があるものの、

結局は偏微分方程式を解いているので、共通のスキームになります。



流体の解析をする場合、おそらく機械流体や熱流体を取り扱う人は1次精度のスキームを使うことは避けると思います。



でも、既存のソフトウェアを使った解析はどうでしょうか?

1次精度スキームを選択できないものはないと思います。



要は数値解が飛んじゃ駄目なんです。

結果が得られないと。

間違ってても、解ければ良いんです。

そんな認識の人は数値シミュレーションを辞めた方が良い。

選択できることは重要ですが、真値からの誤差が大きい、下手したら、解きたい現象を解けていない可能性がある、その認識を持ち続けたいものです。


Common Lisp

実践Common Lispを購入しました。

少ししか読んでませんが、手取り足取り感がハンパないです。

Lisp

Lisper。

カッコいい響きです。
第1言語はLisp。日常からアプリまでLisp。
Emacsしか使わない。

こんな人々はすごいですね。
「自己進化していく言語」を操るウィザードたち。

今日も勉強します。



Visual Studio

Visual Studioが必要になり、会社で買ってもらいました。

いや、正確にはIntel Compiler C/C++が必要だったんですけどね。
統合開発環境がないとなにぶんデバッグとか掛けづらいので。

2012を初めて使ったんですが、インターフェイスがかなり変わってました。
それでもコーディングはしやすいですね。

あと、速い・・・。書き方の問題はもちろんあるのは重々承知ですが、
GCCに比べてワンオーダー速い。
Fortran以上に差がでるようです。
もう、戻れないんだろうな〜。

数値解析のみ考える場合は、コンソールで構わないので、実行速度が速いファイルを吐き出せるコンパイラを使うのが一番手っ取り早いです。

PGIも使ってみたいですね。もちろんMac用ですよ。
なんでCygwin使わなきゃならんのか?





xojo

インターフェイスに慣れないです。。。
でもやっぱりMacでもWindowsでも同じ感覚で使えるのが良いですね。
うまくCSVに吐き出せないのは私の勉強不足です(笑)

なんとか会社に買わせたい。。。

2013年6月8日土曜日

Xojo

移流方程式を修正しつつあります。
今のところ、TVDQUICK、TVDMUSCL。
あとはLeapFrogも予定しておきます。
しばしお待ちを。明日にでも整理します。

技術屋

民主主義、わらりやすい説明、低コスト、精度向上。

低コスト。これは 異論を挟まない。
精度向上。これも異論を挟まない。

ただし民主主義と分りやすい説明。
これは一部同意しかねます。

専門家は、一般市民出来ない近似精度とコスト感覚のバランスを取るのが大事です。
ただ、その説明は非常に難しい。
さぼってる訳ではなく、先人たちが一生かけてきたことを数年で理解し、他人に説明することなんて無理なんです。説明相手も勉強する余裕がないんです。

だから、専門家は必ず自浄作用をもって、大きく間違わない方向の精度を提示した上で、納税者を納得させるのが必要だと思います。

加えて、一般市民(業界外の方々)に対して、しつこく説明すべきだと思います。そうしないと、技術屋は必ず減ります。魅力がない。やったことが、「税金の無駄遣い!」。
それが世論で、政治家がそれをよしとするなら、誰も目指さないし、下手したら、海外に技術が流出します。自国のグランドデザインを国の民間事業をうまく制御できなれば、荒廃がひどいことになります。

それを受容するだけの納税者がどれだけいるのか?
国家のために自分が何が出来るのか?
真剣に考えましょう。

専門家、非専門家、「そしてそれを短絡的に報道するマスコミ」みんな協力しましょう。信号なくなるかもしれませんよ。道路だって平気で壊れっぱなし。「仕方ないよ。お金内んだもん」あり得ないですね。日本という国家がどんな検討ツールとしての武器をもち、一技術者としてどのような判断が必要なのか。関係者の皆さん、お金を引っ張ることだけなく、グランドデザインとして、日本という国家がどうあるべきなのか?考えましょう。その検討のためにどのような勉強が必要なのか、それを勉強するのが専門学校だったり大学だったりします。

みんな得意、不得意は必ずあります。それが最適化されてるとも思いません。

日本のみなさん、選挙に行って、どう考える、悩みましょう。




C/C++

C/C++は書き方がたくさんあります。

Better Cとして。
C with classesとして。
STLをふんだんに使って。

数値解析としては、実行が速くて、可読性に優れ、修正しやすいものが好まれると思います。速度はfortranとC/C++がトップ集団かと思います。

ただし、可読性や修正のしやすさはイマイチ(独断)です。
いや、可読性はFortranであれば、素直に読んで行けば必ずわかるはずです。
時間はめちゃくちゃかかりますが・・・。Cはポインタの概念、C++であればオブジェクト指向の実装方法が分ってないといくら読んでも次に進みません。

解析のコードを0から作るとなると、一旦、ベースの部分をスクリプト言語で書いてから移植するのが良いんでしょうね。

VBAは邪道ですが、結果の確認や出力も簡単なのでオススメですよ。
簡単なGUIも簡単に使えますし。
 その後、高速に実行できる言語に移植します。

 数値解析に依存しない言語にしたいならC/C++でしょう。

なんかソース書く書く言っておきながら晒せてないです。
済みません。

PDE

偏微分方程式はシンプルです。
でも奥が深い。

数値解析では、保存系、非保存系とか出てきます。
数式的には意味は変わりません。

ただし、その結果、保存系、非保存系の基礎方程式如何により、数値解析を得るスキームが変わります。面白い。そして、数値解析を本業としない先生にはその話が通じません。
なぜなら、偏微分方程式を展開しても、数式上等価だからです。

そして、各科学技術分野において、きれいな方程式を用いることは少ないです。
Source Termが必ず入ります。加えて、比例定数が入ると思います。
そのほとんどは陽的な関数で解けません。
摂動法を用いるとうまく表現できることが出来ます。

その比例定数、もしくはプリミティブ変数の関数として表現される変数をうまく近似することが、工学的に意味があると思ってます。各項がどの程度のオーダーなのか、無視し得るのか、場に応じて無視できるのか。
「数値解析」のみに注力されている人はその感覚が分らないと思います。
解析自体は出来ますし、コーディングも非常に速いですが・・・。パラメータ
それは工学的センスに乏しいです。
摂動法は工学的センスにあふれてると思います。

最近でも摂動法を用いた近似解について論文を出している先生がチラホラ見えます。
非常に嬉しくなります。本当に私みたいな人間に対して、こんな考え方があるんだよ!という意識付けになります。摂動法(特異摂動法)を自由に使えると、視野が広がります。
摂動法による解は、電卓叩けば出てきます。これは工学として素晴らしいんです。収束計算しなくていいですし、明確にどの項に注目するのかが明確です。

今後、記述を充実したいと思います。

今日はその1と言うことで、
 ざっくり著者の思いを書いてみました。

XoJoで数値解をグラフ描画を作成したいです。

以上








2013年6月7日金曜日

xojo

Real basic

ポアソンの方程式や格子生成に使えれば良いかなと。


C/C++

数値解析をやる人間(主に科学技術計算)にとって、C/C++のメモリ管理だと思ってます。普通やれよ!とエキスパートに激怒されそうですが、
Fortanは余気にせず、計算アルゴリズムの組み合わせに集中できるところだと思います。

ただし、好きか嫌いかは別。ポインタを理解すれば、あらゆることが出来ると勝手に思ってます。今のところ、第1言語はC/C++、第2言語はPython,
第3言語はRuby。
一番はLISPを理解することだと思います。如何せん「LISPER」です。このように言える人は良い意味で相当マニアックです。

GPGPUならC/C++ですね。組み込みならCになるかと思います。
 0から作るのはあまりないですが、修正を掛ける場合は、Fortran,
C/C++ですね。

新規のウェブアプリなら、PHPとかRailsになるんでしょうね。


流れの解析

流体力学は物理学の範疇では古典です。
コンピュータの演算能力が向上した今では、基礎方程式を(O)記号で支配要因を明確にして、陽的に表現できれば、学ぶ必要もないかもしれません。

加えて、大先輩方はソースコードを書くなんてことをしないと思います。
解析好きの私としては納得がいかないことが多いです。

模型実験でやる検討を数値解析で条件等を詰め、最小の検討ケースで実験をやる。
これこそが数値解析のすごさであり、実験をやることの有意義さだと思ってます。

なにより、お爺ちゃん世代は微分方程式の解を摂動展開なんかをして、陽的に表現してます。これがハンパないです。
絶対に間違いません。微分項の性質を理解しているかです。

私自身はお爺ちゃんになっても、プログラムをひたすらコーディングしてデバッグ作業。
ニュートン物理学の範疇では、解の大まかな性質を持っておくのが大事すね。この感覚こそ工学者の一番の売りだと思います。
zojoのソースコードについては、もうしばしお待ちを。
CSVの吐き出し、グラフ描画がうまく行きません。
 (誰も気にしない)

real basicは楽しいです。楽しいのが一番です。日本でもサードパーティをいっぱい出して欲しいですが。

まったくモッて回し者ではありません。Basicをベートして、多機能で使えます。インターフェイスはプラットフォームで違いが余ありません。
コンソール、デスクトップも同様に使えますし、プラットフォームをシームレス行けるのはすばらしい。

FortanやC/C++には一切かなわないですが、データ処理には絶好調だと思います。

2013年6月5日水曜日

xojo

アクティベーションがどんな形になるのか、イマイチ分ってない。
4月に買ったのだが、xojoに使えるのか?


テキストディタ

みなさんはどのテキストエディタを使うことが多いでしょうか?

私はmiが最初です。
その後、Xcodeに移り、Windowsを使うようになりEmeditor、MIFESも入れてます。

最後に行き着いたのはEmacsでした。Lispカッコ良い。Matzが使ってるからという、かなりミーハーです。Emacsは圧倒的でした。org-mode、howm、非常に使えると思います。

open-junkで簡単にコードの動作チェックできるし。
簡単にshellを呼び出せるのが良いです。

好みのエディタはそれぞれたくさんあって、この中でどれかを否定する気は一切ないです。皆さん、テヅガクをもって作成されていると思いますし、多くがフリー。

この秀丸とかも優秀ですし、Notepad++、サクラエディタもそうです。
ただし、カスタマイズを考えるとVim、Emacsが圧倒的だと思います。学習速度は必要とされますが。。。rubyだとTextMateのようですが、Emacs(xyzzy)でごり押しする予定です。

Lispの記述も増やしたいと思います。

Xojo

出ました。
早速使ってみると、インタフェース周りがかなり変わっています。

びっくりはしませんでした。
こんなものかなと(その時点でライトユーザだとバレる...)。
でも、デバッグのみならフリーのはずです。

Basic系の言語がかなり残っているのは、やはり先人たちが使っていること、
加えてxojoはオブジェクト指向をバリバリ使っているので、かなりコーディングが短くて済みますし学習コストも低いと思います。

オブジェクト指向も学び安いと思います。rubyを勉強したあとにコーディングするとよくわかります。

サンデープログラマーとして、一番良いと思います。
GUIプログラミングでは抜群です。プラットフォーム依存しませんし。

強いて言えば、もう少しリファレンスや書籍を充実させれば、敷居が下がるように思います。Visual Basic(or VBA)を使っていれば入りやすいですし。

予定では線形移流方程式の解法を記述したいと思ってます。
以前の一次風上だけでなく、TVD−Lax-Wendroff法やCIP法、日本人なら外しちゃだめでしょK-Kスキームは予定です。あとはCSVへの吐き出しと、グラフ描画までできれば。

詰まったところについても記述したい意向です。

GUIをネィテブコードで吐き出せて、プラットフォームに依存しないのは非常に素晴らしいと思ってます(回し者ではないです!)。MacもWinもLinuxも使う身としてはありがたいです。これこそサードパーティに良さ。

関係者の皆さんに感謝です。
つたないコーディグになると思いますが、「出来損ないのアマチュア」がバカなりに頑張ってると思っていただければ・・・。

Real Basic5.5の日本語マニュアルが欲しいです・・・。

2013年6月4日火曜日

vim

「i」挿入モード
「v」ビジュアルモード
「esc」ノーマルモード

ビジュアルモードで
「y」:ヤンク
「p」:ペースト

「:w」保存
「:!」シェル
「:q」終了

2013年6月3日月曜日

C

最近、Cを書いてない・・・。
高速の処理にはどうしてもC/C++が必要だから、忘れないようにしないと。

と言いつつLispを入れてます。

最高速はFORTRAN、C/C++。
続いて、Javaとかですかね。
私の独断と偏見に満ちあふれてますが。

スクリプトだとPythonが速いと思います。
rubyの方が若干遅いかなという印象です。

一番カッコいいのはLispだと思います。

今の状況は先人たちの努力に立ってるもの。
皆さんに感謝です。

2013年6月2日日曜日

乱流

NS方程式はそのままだと解けません。
というのは語弊がありますが、現在のスパコンを持ってしても、限られた流れしか解けないはずです。それは流れというものが、非常に小さい時間スケールで変動するためであり、
それを解像するには空間の刻み幅も同様に小さくしなければならないからです。

ただし、その非常に小さい時間スケールの流れの速さ(流速)が問題になることは限られています。もっと大きい時間スケール(数分間の平均とか1時間の平均とか)での速さを知りたい場合がほとんどかと思います。

そんなときにRANSを使います。ある程度のスケールの乱れを解像したいのであれば、LESを使います。

最近、乱れを勉強し直して、楽しくなってきました。
スペクトル解析も使えないといけませんし、摂動法も重要です。

勉強して楽しいと思えるのが良いですね。

Fortran

BLASとか、LAPACKとか、いろいろ入れてみたいですね。
速そうですし。

fortran

カンマの場合は注意が必要。

program main
  integer i
  integer fout
  real(8) a(100),b(100)
  fout = 20
  open(fout,file='test.out')
  do i=1,100
    a(i) = real(i)
    b(i) = real(i)
  end do
  do i =1,100
    write(fout,'(I5,a,e10.3,a,e10.3)') i,",",a(i),",",b(i)
  end do
end program main

2013年5月30日木曜日

Ruby

初学者にとってRubyの良いところ(と勝手に思ってます)。

オブジェクト指向を学ぶのに一番良い(と勝手に思ってる)。
汎用プログラミング言語。
数値解析を出来る。
Emacsを勉強したくなる。
作者が日本人。

もっとたくさんあると思います。
日本発信のプログラミング言語。

使うことで、応援したい。

解析とは

解析解と言うと、私は陽的な関数で表現できるものを想像します。
評価がしやすく、電卓を叩けば解を得られます。
工学的な立場からすると、コストが低く済みます。
ただし、精度の問題があろうかと思います。

「陽的な関数で表現」できる現象は、それほど多くないと思われます。
したがって、ある程度の精度を持つ解を得るために「数値解」が必要になると思います。
離散化して空間積分、時間積分すれば解は得られるでしょう。

それとは別に摂動法があります。
パータベーション。

非常に有用だと思います。そして「オーダー」の概念が重要であり、工学の極みのように思います。

プログラムはものづくりに近いです。
数値解を求めるのは非常に楽しいです。

でも、技術屋としてはちょっと立ち止まり、精度とコストのバランスを考えて、最適な手法を選びたいものです。

摂動法も実は楽しいし、好きだったりします。

「目的は何なのか」
「どんな答えが欲しいか」
「そのためにどの程度の精度が必要か」
「そのためにどれくらいのコストが必要か」
「得られた解に対してどのように評価するか」
「結論は何なのか」
「目的と結論に整合が取れているか」

こんな手順でPDCAを回して、次の世代に明確な課題(目的)を残す。これが技術者としての理想像だと思います。




2013年5月28日火曜日

xyzzy4

やはり、設定とか、使い方とか記した方が、
関係者のモチベーションになるというか、恩返しのような気がしてきました。
Windowsを使ってて、今までで一番すごいソフトだと思いました。
(興奮冷めやらぬ・・・。)

暇を見つけて、記しておこうと思います。

AWKと絡められたら、良いな。

音楽

情熱スパイラルというアルバムを買いました。
ベースがかっこいいです。しかも女性。

ギターのギュイーン(不協和音)も良かったです。
でも、う〜ん、やはりバンドはリズム隊から入ってしまいます。
 ベースの音を拾いながら聞いてしまいます。

チリペッパーズのフリーさんの影響だと思います。
「Give it away」は格好良すぎます。

メロディだとどうしてもピアノの方が好きですね・・・。

2013年5月27日月曜日

xyzzy3

今日設定したのは、
・fortran
・ruby-mode
・outline2

特にoutline2は本当にすごい。関係者に感謝です。emacsユーザなら使いやすいですし、プログラム書くのにも良いです。あんまり期待して無かった分、正直ここまでとは…。
軽い…。落ちない…。

作ってくれた方、メンテされてる方、重ね重ねありがとうございます。

xyzzy2

とりあえずインストール。
後はキーワードとか、wwwのブラウザとか、補完とか出来たら、言うことなし。

2013年5月26日日曜日

xyzzy

普段はNTEmacsを使ってますが、良く落ちます。
Mac、Linuxは落ちないんですが。。。

会社はWindowsなので、xyzzyを使ってみようかと本だけ購入。
でも一番はCommonLispの簡易環境が整うところだと思う。

なるべく、少ないツールを使い込みたいですね。

2013年5月22日水曜日

python

a = []
a.append(2)

言語の使い分け

大学のとき、C言語の単位を落としました。
某大学のコンピュータリテラシーはUNIXなので、WindowsやMacすら触ったことなかった人間が、いきなりターミナルなんかでls、pwd、cd、mkdirなんかで操作を覚えて、muleを使ってました。

数年後、本気でCを勉強しようと思ってポインタの本を買って、すごく世界が広がりました。配列とポインタのつながり、構造体変数のアクセスなんかがわかると面白くなりました。

あくまで科学技術計算の範疇なので、流体の差分法や有限体積法を数値的に解くツールにすぎないので、メモリの管理とかに気を取られたくはないですが、如何せん、ほかの言語に比べて速いです。科学技術計算目的ならFortranとC/C++の2択になると思います。

一方で、出てきた解析結果の整理までFortranやC/C++を使うことはなく、VBAやその他スクリプト言語の利用が多いと思います。FortranやC/C++がヘビー級、スクリプト言語はライト級でしょうか。二つくらいあると良いと思います。
  • ヘビー級(Fortran、C/C++ インテルコンパイラだと激速、PGIは使ってないですが、GPUに計算させるならPGIのようです)
  • ライト級(Perl、Python、Ruby、Pythonのmatplotlibはすごいです)
  • ワンライナー(Awk、ほとんど勉強しなくてもCを知ってればかなり使えるはず)
ただ、科学技術計算に特化すると、Java、PHP、JavaScript等のWeb関連で利用頻度の高い言語には触れないことが多いですね。なんか馬鹿にされている気分なので勉強した上で、敢えて上記のヘビー級とライト級を使いたいです。
AwkのワンライナーはLinux文化が色濃いですが、むちゃくちゃ便利ですね。なぜ広まらないのか不思議です。

遊びだったり、趣味のプログラムならRealBasicを使いたいです。これが仕事に出来たらいいなあと最近思っています。GUIでいろいろ追加できるので便利ですよ。

統合開発環境

例えばRubyを使うとする。
エディタはEmacs、Vimが一番良いと思う。
プラットフォームに依存せずにそれぞれバイナリがありますし。
一度、覚えるとかなり速くなる。

ただ、大規模開発環境になるとエディタだけではかなりきつくなる。
プラットフォーム別に考えると
  • Windows
  • Mac
  • Linux(個人的にはCentOS使ってます)
が大まかに分けられるはずです(FreeBSD等UNIXもありますが・・・。)。どの環境を使うかでツールを使い分けるのはばからしいです。
いろいろ見たが、フリーで考えるとEclipse、Aptana、Netbeansがなるものがあるらしいですが、RubyやPythonだとAptanaが一番良さげの印象です。
WindowsとMacで試しましたが、使用感は余り変わりません。

ところでRealbasicも勉強中ですが、言語および開発環境として個人的に一押しです。
使い捨てのスクリプトも大規模も行けます。vbaをかじってる人は多いと思うので、違和感なく使えます。Webアプリもいけるのでもっと勉強して、この備忘録に追加できれば良いと思ってますよ。

realbasic pure advection

' Sub Actionに記述した線形移流方程式
’f:スカラー量
’fn:スカラー量 (更新後)
’df/dt + u df/dx = 0
'dはパーシャルに読み替えてもらえれば。

 Dim s as string
  Dim i as integer
  dim f() as double
  dim fn() as double
  dim n As  Integer
  dim fout as FolderItem
  Dim txtout as TextOutputStream
 
  s = "Program START!"
  lbResult.AddRow(s)
  Redim f(ixmx)
  Redim fn(ixmx)
 
  For i = 1 to ixmx
    f(i) = 0.0
    fn(i) = 0.0
    If i>10 And i<20 Then
      f(i) = 1.0
      fn(i) =1.0
    end if
    s =str(i) +","+ str(f(i))
    lbResult.AddRow(s)
  next i
 
  For n=1 to 100
    for i =2 to Ixmx
      fn(i) = f(i) - u*dt/dx *(f(i) - f(i-1))
    next i
    for i = 2 to ixmx
      f(i) = fn(i)
    next i
  next n
  For i = 1 to ixmx
    s = str(i) + ","+str(f(i))
    lbResult.AddRow(s)
  next i

2013年5月21日火曜日

スクリプト

# タブ区切りからカンマ区切りに変更
{
  gsub("\t",",",$0);
  print $0;
}

Vim

しばらく、Linuxの勉強もかねて、VIMをやろうと思う。
Ctrlキーを使わないことに違和感が(笑)

Pythonはこっちの方が良いかも。
あとはPHPだな。しばらく開発環境は使わずにシコシコソース書くことにしてみる。

JavaScriptはC++で行き詰まった時に、息抜きになる。
みんな使えるしね。エクセルVBAと同じくらい広まっても良いと思う。

TVD

方程式を直接離散化した方が、分りやすい。
分りやすいとは、余計なスキームの特性を考えなくて良いということだろう。
数値計算の手法を研究する分野ならそれは駄目だけど、ユーザー(数値解析をツールとする人間)からすれば、優先順位は下がる。

その点では、パタンカー先生(香月先生の訳本)の本は、
各微分方程式の各項の解に集中できるので、すばらしいと思う。古いですが。

でもあまりにも数値シミュレーションが当たり前になりすぎて、
本当の基本を分ってる人が少ない。

個人的にはすごくCIP法が好きです。
でもTVDを勉強するとさらに移流方程式に対する理解が深まる上、やっぱりCIPの素晴らしさを改めて実感できる。CIPで数値計算好きになったのもありますが・・・。

こういうことを分った人間が数値解析して、他人に説明しないと、どんどん日本は駄目になると思う。

2013年5月19日日曜日

2013年5月15日水曜日

C/C++

インテルコンパイラが欲しいな。
やっぱり実行が速いしね。

2013年5月14日火曜日

ぼうや

awkの処理がうまくいかない。
なぜだ!

ぼうやだからさ!


2013年5月13日月曜日

ruby

Rubyはオブジェクト指向の勉強にすごく良いと思う。
決して、教育用に限定する訳じゃなく、言語の思想がすごい。
pythonもオブジェクト指向で書けるけど、なんか、違和感が・・・。

おそらく、rubyが「純粋」オブジェクト指向だからだと思ってます。
異論反論たくさんあると思いますが・・・。

でも日本人だからRubyを応援したい。
東京近郊のrubyistのあつまりに参加したいけど、
知り合いが誰もいない・・・。

real basic

たのしくプログラミングしたい。

環境の設定だって、めんどくさいし、思いどおりいかないことは多くある。
そんなところに時間をかけたくない人はたくさんいると思う。
(私はそんなところも含めてプログラムが好きだが)

real basicはフリーじゃないけど、とりあえずインストールしたら、
プログラムをすぐ書けるようになるし、プラットフォームに依存しないし。
WindowsでMacの実行ファイルを作成できるのはすばらしい。

早く、iOS対応してくれたらなあ。

C

Cのコンパイラも欲しいと思っている。
インテル製のやつ。

Cは最近本気で書いてないので、忘れてるんじゃないかと恐怖が・・・。
あ、C++じゃないですよ。純粋なCの方。

ヘッダーは仕様書だ、的なことを聞いて、本当にそうだなと思う。
動的言語はそれがないから、使い捨てに近い使い方しか出来てない。
(いや、ちゃんと書けばいいことだけで動的静的議論をするつもりはないです)



fortran

fortranはgfortranがフリーだけど、
思い切ってIntelのFortranをOS X版を買った。
むちゃくちゃ早いです。

がんがん使って、自分の能力磨かないと。


AWK

プログラミングAWKを買ってみたけど、
すごく面白い。
文法の説明だけでなく、後半のアルゴリズムは思考の訓練になると思う。
名著はやはりすばらしいし、自分の中で長持ちします。

「10日で〜」
「誰でも分る〜」
も否定はしないけど、すぐに飽きてしまう。

Real Basic

Real basic
WindowにPushbuttonを配置。
pusshbuttonをダブルクリックして、Actionの関数で下記ソースを記述

MsgBox "Hello World"

これでボタン押すとHello Worldが出てくる。

GUIプログラミング

GUIのプログラミングは世の中的には普通だろうけど、
科学技術計算ではコンソールが主流だと思う。

real basicはGUIのプログラミングを実施するのに、すごく良いと思う。
VBAをかじっただけだが、何となく書けちゃう。
IDEも優れていると思う。Textでもソース出せるし。

会社で使わせるためのロジックを組み立てなければ。

AWK

AWKのフィルタがすごい。
# test.awk
{
      print $0
}
でファイル登録。シェル上で
gawk -f test.awk input.txt > output.txt
で「input.txt」のデータを読み込んで、「output.txt」にリダイレクトされる。
基本文法がCと類似しているので、使いやすい。