ベクトルの外積を含む微分方程式

home > 掲示板 >

このページのPDF版

サイトマップ

ベクトルの外積を含む微分方程式

yukihiro さんの書込 (2005/06/01(Wed) 09:01)

こんにちはyukihiroです．ボールの運動シミュレーションを行うにあたり，どうにも解けない問題にぶち当たったので，ヒントをいただければと思います．

球が回転しながら空気中を飛んでいるとき，回転ベクトルと球の速度ベクトルとの外積の方向に力を受けます．マグヌス効果というそうです．

初期位置と初速，および回転ベクトルがわかっているとき，任意の時刻における速度と位置を算出したいため，以下の運動方程式を立てました．空気抵抗は速度に比例するものとし，とりあえず風速は0とします．

md<b>V</b>/dt = m<b>G</b> - k<b>V</b> + s<b>V</b>×<b>R</b> ・・・(1) m：質量 <b>V</b>：速度ベクトル t：時刻 <b>G</b>：重力加速度ベクトル k：空気抵抗係数 s：マグヌス効果係数 <b>R</b>：回転ベクトル

これを整理すると d<b>V</b>/dt = a<b>V</b> + <b>V</b>×<b>B</b> + <b>C</b> ・・・(2) a：定数 <b>B</b>,<b>C</b>：定ベクトルと表せますが，一般にこのような微分方程式を解くにはどうしたらよいのでしょうか？

スカラーvの場合は dv/dt = av + biv + c ・・・(3) a,b,c：定数 i：虚数単位

を解くと v = (v ₀ + u)e ^at e ^bit - u ・・・(4) v ₀ ：初速 u：= c / (a + bi) = c(a - bi) / (a ² + b ² ) = (ac - cbi) / (a ² + b ² ) e：自然対数の底

が得られるので，これを(2)式のベクトルの場合にあてはめると e ^bit の乗算は回転とみなせるのではと思い，以下のようにしてみました． <b>V</b> = (<b>V</b> ₀ + <b>U</b>)e ^at M(<b>B</b>t) - <b>U</b> ・・・(5) <b>V</b> ₀ ：初速度ベクトル <b>U</b>：= (a<b>C</b> - <b>C</b>×<b>B</b>) / (a ² + |<b>B</b>| ² ) M(<b>B</b>t)：<b>B</b>/|<b>B</b>|を回転軸とし，|<b>B</b>|tラジアン回転させる変換行列

実際これをプロットしてみると，<b>V</b>と<b>R</b>が垂直な場合しかそれらしい動きを見せませんでした．

四元数をうまく使えばできそうな気がするのですが，答えを導くに至りませんでした．

長々と書いてしまいましたが，要するにお聞きしたいのは(2)の微分方程式を解く方法です．どうか救いの手を・・・

Re: ベクトルの外積を含む微分方程式

崎間さんのレス (2005/06/01(Wed) 17:17)

> 要するにお聞きしたいのは(2)の微分方程式を解く方法です．

外しているかもしれないですが，成分ごとに計算してみてはどうでしょう？

Re: ベクトルの外積を含む微分方程式

Joh さんのレス (2005/06/01(Wed) 17:19)

成分に分けちゃ，ダメでしょうか？

Re: ベクトルの外積を含む微分方程式

yukihiro さんのレス (2005/06/01(Wed) 20:36)

お二方ありがとうございます．

例えば(2)式のもっと単純な場合，a=0，<b>C</b>=零ベクトルの場合を考えると d<b>V</b>/dt = <b>V</b>×<b>B</b> となりますが，成分分解するということは， <b>V</b>のx成分をVxと表すとすると dVx/dt = VyBz - VzBy dVy/dt = VzBx - VxBz dVz/dt = VxBy - VyBx ということでしょうか？

これはどうやって解くのでしょうか？

Re: ベクトルの外積を含む微分方程式

Joh さんのレス (2005/06/01(Wed) 23:12)

積分変数がtですので，とりあえず，右辺をtが見える形に書き直してみてはどうでしょうか．

Re: ベクトルの外積を含む微分方程式

yukihiro さんのレス (2005/06/02(Thu) 08:27)

tが見える形・・・

dVx/dt = Bz∫(dVy/dt)dt - By∫(dVz/dt)dt ということ・・・か・・な・・？？

ごめんなさい．頭が固いもので展望が見えてないです・・・よかったらもう少しだけ噛み砕いてご教示いただければ・・・

Re: ベクトルの外積を含む微分方程式

Joh さんのレス (2005/06/02(Thu) 10:25)

すみません．いい加減なことを書いたかも知れません．Vの関数形は，全然分かってないのですね．

では，行列で表してみるのはどうでしょう？

Re: ベクトルの外積を含む微分方程式

崎間さんのレス (2005/06/02(Thu) 11:13)

yukihiroさん，Johさんこんにちは．横槍すいません．せっかくシミュレーションするのですから，直接数値計算しちゃうのはどうでしょうか．

http://irws.eng.niigata-u.ac.jp/~chem/itou/cemath/cemath08.html

とか

http://www.geocities.jp/supermisosan/rksimultaneousequation.html

が参考になると思います．

Re: ベクトルの外積を含む微分方程式

toorisugari no Hiro さんのレス (2005/06/02(Thu) 19:23)

通りすがりのHiroです．

>これを整理すると >dV/dt = aV + V×B + C ・・・(2) >a：定数 >B,C：定ベクトル >と表せますが， >一般にこのような微分方程式を解くにはどうしたらよいのでしょうか？

両辺を |a| で割って改めて

$t &:= |a|t\\B &:= B/|a|\\C &:= C/|a|$

とし， a<0 で有ることに注意すると Error : LaTeX command contains irregal 2bytes character(s). になります．が定数ベクトルならばの方向に軸をとれば B=(0,0,B_z) とできるので，成分で展開してください．行列 $\hat\Omega$ を

$&1 , -B_z , 0 \\&B_z , 1 , 0 \\&0 , 0 , 1$

で定義するなら $dV/dt = -\hat\Omega V + C$ となり， $U := V - \hat\Omega^{-1} C$ とすれば $dU/dt = -\hat\Omega U$ と簡単になります．成分は自明ですし， x,y 成分は特性方程式（ケリーハミルトン） x^2 + 2x + (1+B_z^2) = 0 より微分方程式 (d/dt)^2 f(t) + 2 (d/dt)f(t) + (1+B_z^2)f(t) = 0 の解となります．

Re: ベクトルの外積を含む微分方程式

yukihiro さんのレス (2005/06/03(Fri) 07:17)

皆さんありがとうございます．

＞Johさん行列で・・・少し考えてみます．

＞崎間さんルンゲクッタ法を検討したこともありましたが，この方法だと逐次的にしか値を算出できないので任意の時刻における速度や位置を定数時間で求めたいというニーズには合わないと思い，使用しない方向でいました．

が，結局アニメーション表示をさせる以上，一定の細かさで段階的に位置を算出する必要があるわけで，それならルンゲクッタ法でもいいかなと思い始めました．これなら微分方程式が複雑になっても (たとえば球の回転速度が空気との摩擦で減衰するとか・・・) あらたにそれを解く処理を追加しなくてもいいわけですし．

あらためて検討させていただきます．ていうか巷の3D物理シミュレーションではこれが一般的なんでしょうかね．

＞Hiroさんなかなか難解ですね^^;

微積分の知識があまりないので，突然特性方程式が出てくるところがよくわかりませんでした．

そこはおいておくとして， x ² + 2x + (1 + B _z ² ) = 0 は虚数解 x = -1±B _z i を持つので，示された微分方程式の解は(大学時代の微積の教科書によると)

f(t) = e ^-t (C ₁ cos(B _z t) + C ₂ sin(B _z t)) C ₁ ，C ₂ ：積分定数

となるわけですが，２つの積分定数はどうやって決定するんでしょう？ t=0において<b>V</b>=<b>V</b> ₀ という１つの初期条件しかないわけですが・・・

Re: ベクトルの外積を含む微分方程式

yukihiro さんのレス (2005/06/03(Fri) 08:29)

追伸座標軸を取り直すというアイデアはいいヒントになった気がします．このアイデア自体はこの「物理のかぎしっぽ」にもありましたね．

「角速度と速度との関係」

http://www12.plala.or.jp/ksp/mechanics/angularVelocityAndVelocity/index.html

これを応用することを思いつかなかった自分の頭の固さを実感します．ていうか少し難解でちゃんと読んでなかったのが悪いのかもしれませんが．

Re: ベクトルの外積を含む微分方程式

yukihiro さんのレス (2005/06/03(Fri) 08:47)

何度もすいません．

> ２つの積分定数はどうやって決定するんでしょう？
> t=0においてV=V0という１つの初期条件しかないわけですが・・・

何となくですが， C ₁ = C ₂ としてかまわない気がしてきました．根拠はちゃんと説明できないですが，何となくそんな匂いがしたもので．

Re: ベクトルの外積を含む微分方程式

おこめさんのレス (2005/06/03(Fri) 11:34)

どうもはじめましてyukihiroさん．理学部の学部3年のものです．記事の引用ありがとうございます(この記事を持ち込んだ者です)．

＞・・・ていうか少し難解

すみません．改訂をしようとは思っているのですが・・・

＞２つの積分定数はどうやって決定するんでしょう？

ですが，これは２つの条件を与えれば良いと思います（これがシュミレーションする内容を示すパラメータの一部では？1自由度あたり2つの積分定数，あと定ベクトル2つと1つのスカラーですから自由度3の運動をシミュレーションすれば計13個のパラメータでしょうか？）．他のパラメータを固定したとすると，1自由度あたり２つの独立な運動条件のみが，運動のバリエーションを与えていていると思います．仮にこれが任意の状況で完全に決まってしまったら，どんな状況でも同じ運動をしてしまうという何ともおかしな状況になってしまいそうです．

たとえば，初期速度と初期位置のみで運動（軌道）が完全に決まると思えば良いのではないですか（初期位置，初期速度が違うのに同じ軌道を描くのが正しくないのはシミュレーションせずとも分かると思います）？

Re: ベクトルの外積を含む微分方程式

toorisugari no Hiro さんのレス (2005/06/03(Fri) 19:15)

> 突然特性方程式が出てくるところがよくわかりませんでしたえっと，たとえば

$(d/dt)x &= a x + b y\\(d/dt)y &= c x + d y$

という微分方程式を考えます．を消去してだけの式にしたいのですが，結構大変です．でも行列 $\hat A$ を

$a ,b \\c ,d$

とすれば上の式は

$(d/dt)(x,y)^t &= \hat A (x,y)^t$

とかけます．( ~^t は転置) ここで，特性方程式は $\hat A^2 -(a+d)\hat A +(ad-bc)\hat 1 = 0$ ですから，

$\hat A^2 (x,y)^t-(a+d)\hat A (x,y)^t + (ad-bc) (x,y)^t&= 0$

であり，結果

(d/dt)(x,y)^t - (a+d)(d/dt) (x,y)^t + (ad-bc) (x,y)^t&= 0

がえられます．

あと，初期条件ですが，で2個，で2個の計4個でますが，方程式に代入すれば2個減るので最終的には2個となり， V_x^0 , V_y^0 を与えれば完全に決まります．

Re: ベクトルの外積を含む微分方程式

おこめさんのレス (2005/06/03(Fri) 22:12)

上の解説は特性方程式自体の意味が分からなければ厳しそうです（よく読んでませんが，はじめから連立方程式を考えているみたいなので）．簡単のために定数係数線形微分方程式

$y^{\prime\prime}(x)+ay^{\prime}(x)+by(x)=c(x)$

の解を求めるときを考えます．簡単のためにまず c(x)=0 の場合を考えて，解を初等関数の一般形 $y=e^{kx}x^{\omega x}$ と仮定します（線形性から考えて妥当な仮定だと思います）．これを上の方程式(ただし C(x)=0 )に代入すれば，特性方程式が得られると思います．数学的に正確な説明かは分かりませんが，とりあえず正しい結果は得られます．あとはこの同次形の解から定数変化法で求めるなりすれば上の方程式の一般解が求まると思います．

的外れなら書き込みならば，どうぞ無視してください．お願いします．

Re: ベクトルの外積を含む微分方程式

Joh さんのレス (2005/06/04(Sat) 00:43)

線形なので，微分は，演算子の掛け算だと思って良いのです．そう思えば，ただの二次方程式と，似たようなものなわけです．行列で，ケーリー＝ハミルトンの公式というのがあったのは大丈夫でしょうか．

Re: ベクトルの外積を含む微分方程式

yukihiro さんのレス (2005/06/04(Sat) 18:15)

皆さんありがとうございます．

特性方程式自体は多分高校レベルの話でしたね．すっかり失念してました．それと連立微分方程式の関係もおかげさまで理解できました．

初期条件についても，Hiroさんのおっしゃるとおり連立方程式に代入することにより，４つの積分定数を確定することができました．

おかげさまで，期待通りの動作をさせることが・・・できました，と言いたいところですが，初期条件によって微妙に変な動きをするので現在デバッグ中です．

上記積分定数を求める過程でどっか計算違いをしたかな？などと思ってますが，とりあえず都合により今日明日は進展しそうにないので，また明後日あたりに報告したいと思います．

Re: ベクトルの外積を含む微分方程式

yukihiro さんのレス (2005/06/07(Tue) 15:18)

遅くなりましたが，結果報告させていただきます．

計算が多少複雑なので・机上の計算にミスがないか・それを正しくコードに落とせているかというチェックに手間取り，デバッグにかなり苦労しましたがようやく期待通りの動作をさせるに至りました．苦労した分喜びも一入といったところでしょうか^^;

ちなみにバグの内容は，やはり積分定数の計算違いもありましたが，それに加えてBの方向をZ軸にとるという座標変換についてもミスってましたorz

ここで結論だけまとめておきたいと思います．

> d<b>V</b>/dt = a<b>V</b> + <b>V</b>×<b>B</b> + <b>C</b> ・・・(2)
> a：定数
> <b>B</b>,<b>C</b>：定ベクトル

の解は以下のようになります．

以降すべて<b>B</b>の向きがZ軸正方向になるよう座標を取り直したものの結果で，<b>B</b>≠零ベクトルとします．またx,y,zはベクトルのx,y,z成分を表します．

Vz(t) = (V ₀ z + Cz / a)e ^at - Cz / a (a ≠ 0) Vz(t) = tCz + V ₀ z (a = 0)

Vx(t) = (D ₁ cosbt + D ₂ sinbt)e ^at - Ex Vy(t) = (D ₂ cosbt - D ₁ sinbt)e ^at - Ey b = |<b>B</b>| D ₁ = V ₀ x + Ex D ₂ = V ₀ y + Ey <b>E</b> = <b>C</b>Ω ^-1 Ω = a-b0 b a0 0 0a ※ここで示すΩはHiroさんの示されたΩをa倍して転置したものになります．

V(t)を積分した位置ベクトルP(t)については以下のようになります．

Pz(t) = (V ₀ z + Cz / a)e ^at / a - tCz / a - (V ₀ z + Cz / a) / a + P ₀ z (a ≠ 0) Pz(t) = t ² Cz / 2 + tV ₀ z + P ₀ z (a = 0)

Px(t) = {(aD ₁ - bD ₂ )cosbt + (bD ₁ + aD ₂ )sinbt}e ^at / (a ² + b ² ) - tEx - (aD ₁ - bD ₂ ) / (a ² + b ² ) + P ₀ x

Py(t) = {(aD ₂ + bD ₁ )cosbt + (bD ₂ - aD ₁ )sinbt}e ^at / (a ² + b ² ) - tEy - (aD ₂ + bD ₁ ) / (a ² + b ² ) + P ₀ y

最後に座標軸を元に戻して完了です．

ある程度予想してましたが，回転速度(orマグヌス効果係数)を上げると非常に面白い軌跡を描きますね．

崎間さん，Johさん，Hiroさん，おこめさん，本当にありがとうございました．

Re: ベクトルの外積を含む微分方程式

おこめさんのレス (2005/06/07(Tue) 19:20)

計算は追っていませんが，うまくいったようで良かったです．

No.5448を読んで，僕は定数を完全に決めようとしているのかと誤解してしまいました．それで僕は，勘違いな事を2回ほど書いていただけです．どうもご迷惑をおかけしました．

Re: ベクトルの外積を含む微分方程式

やかんさんのレス (2005/06/08(Wed) 15:07)

yukihoroさん，はじめまして．解決したようでよかったですね．ベクトル，行列の演算は私にはお手上げなのですが最近読んだ岩波科学ライブラリーの”魔球を作る”という書籍に野球の変化球の原因は，ボールの回転による，流体である空気との間に生じるマグナス力だと書いてあったので，きっとその力のことかなあ，なんて想像しました．その本には”直球は変化球？！”と書いてあって直球はバックスピンの回転がかかっているため上向きにマグナス力が働き空気中での放物線（て言わないか・・）よりも浮いている（落ちていない）そうですね．一方フォークボールは回転がほとんどなく，自然の落下にちかいから，フォークボールは変化球ではないそうです．混ぜて投げると，フォークボールの方が ”落ちる球”とバッターには見えるとか・・・．ところでyukihiroさんがお書きの指数や添字はどうやっていらっしゃるのでしょうか．texではないような・・・．

Re: ベクトルの外積を含む微分方程式

yukihiro さんのレス (2005/06/08(Wed) 17:17)

＞おこめさん迷惑なんてとんでもないです．とても勉強になりました． No.5448は気にしないでください．私のいい加減な嗅覚が露呈した恥ずかしいレスなので;;

＞やかんさん

> 最近読んだ岩波科学ライブラリーの”魔球を作る”という書籍に

これですね．

http://www.athome.co.jp/academy/physics/phy05.html

著者のインタビュー記事です．

私もいろいろ調べているときにこのページにも辿り着き，興味深く読ませてもらいました．あんなちょっと出っ張ってるだけのシームがそんなに軌道に影響を与えるものなんですね．巨人の星の消える魔球も可能などと豪語されてますが，信じていいんですかねｗ書籍のほうもちゃんと読んでみたいところです．

> ところでyukihiroさんがお書きの
> 指数や添字はどうやっていらっしゃるのでしょうか．texではないような・・・．

「利用上の注意」に書いてありますよ^^;

http://hooktail.maxwell.jp/cgi-bin/yybbs/yybbs.cgi?room=room1&mode=howto

私の計算式，読み辛くて申し訳ないです．本当はTeX数式を使いたかったのですが，未経験者には少し敷居が高くて・・・この掲示板，プレビュー機能をつけていただくわけにはいかないかな，とかわがままを言ってみたり．

Re: ベクトルの外積を含む微分方程式

やかんさんのレス (2005/06/08(Wed) 19:56)

＞これですね．

http://www.athome.co.jp/academy/physics/phy05.html

著者のインタビュー記事です．

うわー，この先生，結構本の中身しゃべっちゃってますね．きっと学術一筋で，儲けがどうこうなんて考えないのかな :)

＞私もいろいろ調べているときにこのページにも辿り着き，興味深く読ませてもらいました．

すみません・・，ご存知とは知らずに，私ちょっと知ったかぶりしちゃったみたいです．

＞巨人の星の消える魔球も可能などと豪語されてますが，信じていいんですかねｗ

そっかー，これがあるから購買欲がまだまだそそられますね (^^;

＞書籍のほうもちゃんと読んでみたいところです．私は理解できなかったのですが，結構数式とか多いので，お読みになられるときっと面白いのではないかと思います．

＞「利用上の注意」に書いてありますよ^^;

http://hooktail.maxwell.jp/cgi-bin/yybbs/yybbs.cgi?room=room1&mode

=howto 本当はTeX数式を使いたかったのですが，未経験者には少し敷居が高くて・・・

あっ，そうだ！tex以前に使えるようにしていただいたんだった・・．私としたことが．でも簡単な表記だったらなんと私でもできるのだから，yukihiroさんだったら絶対使えちゃいますよ． :)

[home] [掲示板] [ページの先頭]