Previous | ToC | Next |
赤木: というわけで、プログラムとかの解説ね。この章が終わればあ と一つだからさくさくいきましょ。
学生S: はあ、、、(プログラム書いたのも中身の説明するのも僕なの に何がさくさくだよこのやろう、とか思ってるかどうかは内緒)
赤木: とりあえずリスト出して。
学生S: これですね:
# february.rb require 'vector.rb' include Math
def fabs(x) x>0 ? x : -x end def cl(alpha) stallangle=15 maxangle=45 clcoef = 0.05 alpha = atan2(sin(alpha),cos(alpha)) x = alpha*180/PI xabs = fabs(x) if xabs <= stallangle y=xabs*clcoef elsif xabs < maxangle y=stallangle*clcoef*0.4 else y=stallangle*clcoef*0.4*(90-x)/(90-maxangle) end x> 0 ? y : -y end
def cd(alpha) alpha = fabs(alpha) % PI alpha = PI-alpha if alpha > PI/2 x = alpha*180/PI 0.01+x*x*0.0005 end def cds(alpha) alpha = fabs(alpha) % PI alpha = PI-alpha if alpha > PI/2 x = alpha*180/PI 0.01+x*x*0.0005 end def alphaexact(beta) alpha=-beta 100.times{ x = alpha+beta direction=atan2(cl(x),cd(x)) alpha -= (direction-alpha)*0.003 print( alpha, direction, "\n")
} alpha end class Airplane
attr_accessor :m, :sw, :ss, :lw, :ls, :moi, :x, :v, :omega, :psi, :beta, :gamma def theta atan2(@v[1],@v[0]) end
def alphaw @psi - theta end def alphas @beta + alphaw - @gamma*cl(alphaw)+@omega*@ls/@v.abs end
def fwl 0.5* $rho*v.sqr*@sw*cl(alphaw) end def fwd 0.5* $rho*v.sqr*@sw*cd(alphaw) end def fsl 0.5* $rho*v.sqr*@ss*cl(alphas) end def fsd 0.5* $rho*v.sqr*@ss*cds(alphas) end def ev @v*(1.0/@v.abs) end
def nv e=ev [-ev[1],ev[0]].to_v end def fw nv*fwl-ev*fwd end
def fs nv*fsl-ev*fsd end def fg [0,-1].to_v* $g*@m end def vdot (fg+fw+fs)*(1.0/@m) end
def lwvector [cos(@psi),sin(@psi)].to_v*@lw end def lsvector [-cos(@psi),-sin(@psi)].to_v*@ls end def omegadot ((lwvector^fw) + (lsvector^fs))*(1.0/@moi) end end $rho=1 $g=10
def solve(dt,tend,glider) t = 0 ta=[] xa=[] ya=[] vxa=[] vya=[] psia=[] alphaa=[] omegaa=[] while t < tend gtmp = glider.clone h = dt*0.5 gtmp.x += glider.v*h gtmp.v += glider.vdot*h gtmp.psi += glider.omega*h gtmp.omega += glider.omegadot*h
glider.x += gtmp.v*dt glider.v += gtmp.vdot*dt glider.psi += gtmp.omega*dt glider.omega += gtmp.omegadot*dt t += dt
# print t, " ", glider.x.join(" ")," ", glider.v.join(" ")," ", # glider.psi," ", glider.omega,"\n" ta.push(t) xa.push(glider.x[0]) ya.push(glider.x[1]) vxa.push(glider.v[0]) vya.push(glider.v[1]) psia.push(glider.psi) omegaa.push(glider.omega) alphaa.push(glider.alphaw) end $result=[ta,xa,ya,vxa,vya,psia,omegaa,alphaa] end
glider = Airplane.new glider.m=0.02 glider.sw = 0.04 glider.ss = glider.sw*0.1 glider.lw = 0.0 glider.ls = 0.3 glider.moi = 2e-4 glider.beta= -0.05 glider.gamma= 0.1 glider.x = [0,0].to_v glider.v = [7,0].to_v glider.psi = 0 glider.omega=0
default_glider= glider.clone unless File.exist?("february-fig2.ps") glider=default_glider.clone printer february-fig2.ps/vcps square viewport 0.1 1.1 0.8 1 limit 0 70 -6 0 expand 1.5 box relocate 35 -8.5 label X relocate -9 -3 label Y solve(1.0/256, 10,glider) xdata $result[1] ydata $result[2] conn
viewport 0 1 -1.5 -0.5 limit 0 10 -0.1 0.1 box relocate 5 -0.18 label T relocate -1.8 0 label \\gq,\\ga xdata $result[0] ydata $result[7] conn ydata $result[5] conn pgend psfix end
unless File.exist?("february-fig2b.ps") glider=default_glider.clone printer february-fig2b.ps/vcps def cd(alpha) 0 end def cds(alpha) 0 end
square viewport 0.1 1.1 0.8 1 limit 0 70 -2 2 expand 1.5 box relocate 35 -4 label X relocate -9 1 label Y solve(1.0/256, 10,glider) xdata $result[1] ydata $result[2] conn viewport 0 1 -1.5 -0.5 limit 0 10 -0.05 0.15 box relocate 5 -0.13 label T relocate -1.9 0.1 label \\gq,\\ga xdata $result[0] ydata $result[7] conn ydata $result[5] conn
pgend psfix end
赤木: で、どこから話を始める?
学生S: クラスの説明とかしてもあれですよね?とりあえずグラフ書く ところからにします。まず
glider = Airplane.new glider.m=0.02 glider.sw = 0.04 glider.ss = glider.sw*0.1 glider.lw = 0.0 glider.ls = 0.3 glider.moi = 2e-4 glider.beta= -0.05 glider.gamma= 0.1 glider.x = [0,0].to_v glider.v = [7,0].to_v glider.psi = 0 glider.omega=0
default_glider= glider.clone
ここまでは初期設定です。変数名は基本的には本文の表にまとめたのと同じで、 それぞれ Airplane クラスの インスタンス変数で、外からアクセスできるよ うにしたものです。で、これを使い回すのでコピーを default_glider にとっ ておきます。以下が本体です。
unless File.exist?("february-fig2.ps") glider=default_glider.clone printer february-fig2.ps/vcps square viewport 0.1 1.1 0.8 1 limit 0 70 -6 0 expand 1.5 box relocate 35 -8.5 label X relocate -9 -3 label Y solve(1.0/256, 10,glider) xdata $result[1] ydata $result[2] conn viewport 0 1 -1.5 -0.5 limit 0 10 -0.1 0.1 box relocate 5 -0.18 label T relocate -1.8 0 label \\gq,\\ga xdata $result[0] ydata $result[7] conn ydata $result[5] conn
pgend psfix end
solve 以外は全部グラフの枠とかプロットとかしてるだけなので、 solve の ほうをみてみます。
def solve(dt,tend,glider) t = 0 ta=[] xa=[] ya=[] vxa=[] vya=[] psia=[] alphaa=[] omegaa=[] while t < tend gtmp = glider.clone h = dt*0.5 gtmp.x += glider.v*h gtmp.v += glider.vdot*h gtmp.psi += glider.omega*h gtmp.omega += glider.omegadot*h glider.x += gtmp.v*dt glider.v += gtmp.vdot*dt glider.psi += gtmp.omega*dt glider.omega += gtmp.omegadot*dt t += dt ta.push(t) xa.push(glider.x[0]) ya.push(glider.x[1]) vxa.push(glider.v[0]) vya.push(glider.v[1]) psia.push(glider.psi) omegaa.push(glider.omega) alphaa.push(glider.alphaw) end $result=[ta,xa,ya,vxa,vya,psia,omegaa,alphaa] end
このプログラムでは、微分方程式を解きながら線を引くのではなくて、解いた 答を配列に全部しまってそれを呼び出し側に返す、というふうにしてみました。 何故そうしたかは良く憶えてないですが、多分、線で引こうとするとプログラ ムがちょっと面倒になると思ったからです。線を引くには前の値をとっておい て、 2 点を結ぶ必要があるのですが、パネル何枚も使うとかするとその辺が ややこしいですから。
データファイルに落としてもいいのですが、そんなに計算時間がかかるわけで もないし、データの大きさもしれてるので配列にそのままいれました。配列に しまってるのは
xa.push(glider.x[0])
とかのところです。数値積分は2次のルンゲクッタで
の形の公式を使いました。
ではなくて
の値を gtmp
という Airplane クラスのインスタンスにしまうようにしてます。
赤木: なるほど。 2 次にしたのはなんか理由があるの?
学生S: 1次だとさすがに精度がアレだし、4次とかは書くのがちょっと面倒になるので。
赤木: まあ、そんなに長い積分するわけでもないから、これくらいで 十分よね。微分方程式の右辺の計算は?
学生S: これはまあ式の通りなので、説明することもないです。
赤木: んなこといっても読者にはわかんないわよ。 vdot は?
学生S: えーと
def vdot (fg+fw+fs)*(1.0/@m) end
ですね。 fg は重力で
def fg [0,-1].to_v* $g*@m end
だから、下向きの
というだけです。 fw は
def fw nv*fwl-ev*fwd end
で、 nv, fwl は
def nv e=ev [-ev[1],ev[0]].to_v end
あれ、 e って使ってないですね。 ev は速度方向の単位ベクトルで
def ev @v*(1.0/@v.abs) end
fwl に戻ると
def fwl 0.5* $rho*v.sqr*@sw*cl(alphaw) end
で、
で、揚力です。抗力のほうも同じようなものです。
赤木: 関数が大体式1つに対応してるのね。まあ、わかりやすいかも。 でも、なんか無駄が多くない? cl とか cd とか何回も何回も計算されるし。
学生S: 無駄は多いと思います。とりあえず正しく動くプログラムを早 く作るほうを優先にして速度はまあ、、、というつもりで作ったので。色々す れば 10 倍くらいは速くできるかもしれませんが、まあ、見てるうちに計算終 わるので。
赤木: んー、そうね、なんかもったいない気もするけど、こういうのが <web>pitecan.com/articles/Bit/Fugo/fugo.html|富豪的プログラミン グ</a>というのかしら。
学生S: えー、そういう言葉があるんですか?
赤木: あるわよ。これ作者の高校の先輩が作った言葉なの。まあ、今 回はこんなところかしらね。
Previous | ToC | Next |