積算
 | ||||||
| 速度の積算 | ||||||
 | ||||||
| 速度を積算したイメージ図 |
| A | B | C | D | E | F | |
| 1 | No | 速度 | 秒数 | 区間距離 | 距離累計 | 時間累計 |
| 2 | 1 | 速度１ | 時間１ | =B2*C2 | =D2 | =C2 |
| 3 | =A2+1 | 速度２ | 時間２ | ↓ | =E2+D3 | =F2+C3 |
| 4 | ↓ | ： | ： | ↓ | ↓ | ↓ |
Ｂ列Ｃ列には「速度○○で○○秒移動」の実際のデータを順番に入れます。
Ｄ列は、その「○○秒間の移動距離」です(単に速度＊時間)。
これをＥ列で累計していきます。もちろん、最終的なトータルだけでよければ、=Sum()でもＯＫです。
表では「○○秒」のところがすべて２秒です。
この作表ならもちろん５秒などでもＯＫですが、実際に必要なるデータ処理では「○○秒おきに計測したデータ」という場合が一般的です。
最初から「○○秒おき」がわかっていれば(それが普通)、わざわざ秒の欄をつくったりしなくても十分計算できます。
ちなみに、この計算は当たり前なようで、計測や制御の分野などでよく使われています。
数学的には積分に近いため、後述のような数値演算に用いられます。
さて、実際にこの表のように、速度がすぱっと変わるかというと物理的にはあり得ません。
おそらくは、右図の破線のような感じで速度が徐々に増え、一定になり、徐々に減る、というような加減速の過程を２秒おきにはかったら、◆印のデータが得られた、ということです。
しかし、データはこれしかないので、とりあえず、水色の階段状の速度変化を仮定して、縦長の短冊にばらして加えていったのが上の表です。
どうみてもずれがありますが、細かいことは後述します。
蛇足：
工学部だったら「秒速○○メートル」といったときに、「あ、あのくらいか」と思い浮かぶのが望ましいです。
ちなみに、歩く速度が 1[m/s]強(時速4キロ)です。世界最速の人が走って10[m/s](時速40キロ弱)。
差分
 |
| 差分 |
| A | B | C | D | E | |
| 1 | 時刻 | 位置 | 移動量 | 経過時間 | 速度 |
| 2 | 時刻１ | 位置１ | |||
| 3 | 時刻２ | 位置２ | =B3-B2 | =A3-A2 | =C3/D3 |
| 4 | ： | ： | ↓ | ↓ | ↓ |
時刻とともに位置のデータが与えられているとき、単に、位置と時刻の差をとって、割るだけです。
これも、通常は「一定の時間間隔で計測」することが多く、その場合は経過時間が一定となるため、割り算も定数ですみます。
この計算もまた、計測制御の分野で実際によく使われる演算で、微分の代わりにも用いられます。
なお、差や速度の結果である、Ｃ〜Ｅ列で、２行目があけてあります。
これは趣味の問題で、３行目以降を上に詰めてもかまいません。
こういう計算をするときには、「時間が下に向かって進む」ため、「左斜め下＝気分的に時間が後ろ＝未来」の値を使うような印象があるため、このような形にしています。
数値積分
先ほどの積算の計算をする際、積算すべき対象を自分で関数で用意することで、その関数の定積分(に近い値)を計算することができます。
この手法を数値積分、積分の数値(計算)解と呼びます。
一方、微分積分学の講義で習うような、数式のやりくりだけで、最後に値を入れれば結果が出るような結果を解析解といいます。
基本的に解析解は数式変形上の正しさを保って行うため、値も正確です。
しかし、数値解は普通は後述のように誤差が生じます。
誤差は出ますが、解析解のように悩まなくても、そもそも解析解がない場合でも、そこそこ答えが得られることが最大の利点です。
実際に、ここでは積分してみましょう。
 |
| y=2xの数値積分(0.1step) |
 |
| y=2xの数値積分(0.1step) |
 |
| y=2xの数値積分(0.01step) |
を積分していってみます。
答えはいうまでもなく、
です(積分定数＝０)。
ここでは、先の例での時間にあたるｘを一定の幅で増加させ、そのときのｙを数式で与えます。
それ以外はいっしょです。
| A | B | C | D | E | F | |
| 1 | ||||||
| 2 | Δx | 0.1 | ||||
| 3 | ||||||
| 4 | No | x | y=2x | yΔx | 積算 | y=x*x |
| 5 | 0 | =A5*$C$2 | =B5*2 | =C5*$C$2 | =D5 | B5*B5 |
| 6 | =A5+1 | ↓ | ↓ | ↓ | =E5+D6 | ↓ |
| 7 | ↓ | ↓ | ↓ | ↓ | ↓ | ↓ |
この作表でのポイントは以下の通りです。
- ｘの増加分Δｘを特定のセルに置き、絶対参照($C$2)で使うようにする。
→Δｘをあれこれ変えて結果を確認できる。 - ｘもΔｘをもとに作るようにする。B5=0として下向きにΔｘを加えていくのもあり。ここではNoをもとに作っている。
- 積分したいｙをｘ(B列)の数式で作る。
- ｙΔｘをつくる。
- Ｅ列で積算している。計算上はE6=E5+C6*$C$2のように、Ｄ列なしで計算しても全く問題ない。
- 補足：Ｆ列は比較のために、本来あるべき(解析値)ｙ＝ｘ＊ｘをおいている。
- 補足：グラフは積算値Ｅ列と解析値Ｆ列を表示。
- 補足：１行目、３行目の空白はなんとなく見やすくするためのものであり(Δｘが表のパラメータであることを強調)、計算上は意味がない。
このずれが、数値積分をする際に避けられない誤差です。
単純には、Δｘを小さくすれば、誤差は減ります(右図３番目)。
Δｘ＝0.1のときはｘ＝3.0までの積分に当たるx=2.9(※)のところで8.70、これが、Δｘ＝0.0.1のときは8.97と、誤差が0.3から0.03になります。
単純な傾向としては、Δｘを１０分の１に細かくすると誤差は概ね１０分の１になります。
※ｘ＝0.0からｘ＝3.0までを３０等分した状態なので、３０個目のｘ＝2.9(x=2.9〜3.0の区間)がx=3.0までの計算値になる。
 |
| 数値積分の誤差 |
そもそも積分のおおざっぱな解釈は、「関数ｙ＝ｆ(ｘ)と、ｘ軸に囲まれた範囲の面積」です(y＞0ならそのもの)。
右図は「解析値＝曲線の下側の面積＝ピンク」と「数値積分＝短冊の面積の合計＝水色」を図にしたものです。
この図では、曲線が増加しているときには、水色はピンクに足りていません。つまり、少なめに結果が出てしまいます。
逆に、減少時には水色がはみ出しているので、合計値は多めになります。
アップダウンの激しい関数の場合は、適当に帳消しになることが期待されますが、単調増加の場合には不足します。
ここで、短冊の横幅を細くしてみます。明らかにピンクの不足領域も、水色のはみ出しも小さくなっています。
これがΔｘを小さくしたときに誤差が減る理由です。
 |
| 数値積分の誤差と計算方法 |
右図の一番上は「短冊の左端のｘで関数値を使う」です。この場合、先ほど同様、増加時に不足、減少時に過剰となります。
２番目の図は「短冊の右端のｘで関数値を使う」で、逆に、増加時に過剰、減少時に不足します。のでこの二つに大差はありません
３番目の図は「短冊の中点の値を使ってみた」です。なんとなくマシになる感じですが、頂点のところで誤差が出やすくなるなど、抜本的には改善しません。
４番目の図は「短冊をやめて台形にした」です。これは明らかに効果がありそうです。
実際にその効果が認められて「台形積分」と呼ばれる計算方法です。
計算量が増えそうですが、「１本目の短冊の右の高さ(y)」と「２本目の短冊の左の高さ」は同じ値です。ので、
積分＝（１左＋１右）Δｘ/2＋（２左＋２右）Δｘ/2＋..（ｎ左＋ｎ右）Δｘ/2
＝（Δｘ）｛１左/2＋１右(=２左)＋２右＋...(n-1)右＋ｎ右/2｝
となります。もとが、
積分＝（Δｘ）｛１左＋２左＋...＋ｎ左｝
なので、両端がすこし変わるだけです。途中の値がいらないなら、計算量がほとんど変わらないのに精度がよくなるという便利な方法です。
難点は、両端の例外部分(/2)を作り間違う可能性がある、というところです。
いまどきはコンピュータの速度が劇的に速いので、力任せにΔｘを細かくしてすませる、も一つの手です。
 |
| xsinx 0〜π の積分 |
です。やれば、高校クラスの数学知識でも答えはでますが、ぱっとみて「積分したくねー」と思う程度の積分です。
まず、解析解を求めて置きます。
![\int_0^\pi x\sin x dx&=&\left[x(-\cos x)\right]_0^\pi-\int_0^\pi(1)(-\cos x)dx\nonumber\\&=&(\pi-0)+\int_0^\pi\cos xdx\nonumber\\&=&\pi-[\sin x]_0^\pi=\pi-(0-0)=\pi](eqn/spreadsheet2.bhtml.eqn24.gif)
経過の理解は置いておいて、結果はπ=3.14159..です。
次にこれを数値計算できるように変換します。
積分は、先の積算の計算で、横軸の幅を限りなく細かくしたものに相当します。
限りなく狭い横幅dxを関数の値(xsinx)に掛けて積んでいったものです。
計算上は限りなく狭く、という訳にはいかないので、ここでは０〜πを１００等分してみます。
| A | B | C | D | |
| 1 | ||||
| 2 | Δx | =Pi()/100 | ||
| 3 | ||||
| 4 | No | x | y=xsin(x) | 積算 |
| 5 | 0 | |||
| 6 | 0 | =A6*$C$2 | =B6*Sin(B6) | =D5+C6*$C$2 |
| 7 | =A6+1 | ↓ | ↓ | ↓ |
| 8 | ↓ | ↓ | ↓ | ↓ |
積算結果は3.1413。解析解では3.1416なので、まあまああってます。
この分割を細かくしていけば、当然精度はよくなりますが、表計算の上で手頃にできるのは１００くらい、がんばっても実用になるのは１０００分割くらいです。
それ以上細かくして計算精度を上げるには、計算法を見直す(台形積分をつかってみるとか)ことも一つですが、プログラムを書けば、このくらいの積分なら100万分割もわけなく計算できます。これはプログラミングのところで今後やることになります。
この計算は先に述べたように、解析解を求めるのが面倒/不可能(自分には、を含む)なときにそれらしい値を出してくれます。
加えて、「解析解を出せたけど、正しいかが疑わしい」という場合に、検算として使えますし、予めこのあたりの値と目安をつけておくと解析しやすい場合もあります。
そういう「セカンドオピニオン」としての数値解、覚えておくと便利です。
すくなくとも、私の場合は、大学時代には数学の宿題や、研究で必要な数式の確認のため、よく、この手を使っていました。
確実に、着実に結果を求めていく必要がある場合、「常に２系統以上の答えの出し方を考え、結果の一致を確認する」ことは非常に重要です。
余裕があればぜひやってみましょう
- 分割１０００にしてみる。
→$C$2を/1000書き換え、数式を１０００行くらいまで増やすだけ - y=xcosxを積分してみる。
→C列の数式をさしかえ。解析解は２？ 
を−１〜１で積分してみる。
xを−１〜１になるように＆ｙ書き換え。√aはSqrt(a)で求まる(=Sqrt(1*1-B?*B?))。
結果はπ/2=1.57くらいになるはず。これ、何の積分？＜x,yでグラフを書けば一発でわかる？。- 台形積分をやってみる。
→多少計算部分を書き換え
シミュレーション
シミュレーションとは
機械におけるコンピュータの応用事例で多いものにコンピュータシミュレーションがあります。
シミュレーションとは、検討を加えたい対象の挙動を決める各種方程式を元に、順次計算を行っていき、その結果を評価するものです。
これにより、実物での実験をすることなく、様々な結果を得ることができるほか、実物では計測し得ないデータを得ることもできます。
たとえば、
- 流体シミュレーション：自動車や飛行機の周りの空気の流れなど
- 機械的シミュレーション：材料の変形、振動など。加工の模擬なども
- 機構シミュレーション：メカの動きを計算、部品の衝突などが起きないかどうか
- ロボット制御シミュレーション：ロボットを模擬し、それを制御するにはどうするかを検討
- 回路シミュレーション：設計した電子回路が適切に動作するかを検証
- 対象を的確に表す数式
- 動作の条件(多くの場合、計算する領域の端をどうするかが重要で、境界条件と呼ばれる)
- 計算の細かさ
個人的な感想からすると、流体シミュレーションはかなり現実に近い結果が出るようですが、ロボット系、特に歩行ロボットは地面との断続的な接触部分をどうするか、滑りをどうするかが完璧には計算できず、往々にして実物と違う結果がでるようです。
ここでは、シミュレーションを実例で確認してみます。
空気中の自由落下
物体の運動は、
という運動方程式で表されます。言葉で書き直すと、
物体にかかる力＝質量×加速度
です。加速度は、速度の瞬間的な時間変化、つまり位置を時間で２回微分したものになります。
逆にみると、
加速度は力を質量で割ったもの、となります。
ものの運動のシミュレーションの基本はこの式になります。
時々刻々と物体にかかる力(重力とか駆動力とか)を求め(or入力し)、加速度にします。
これを１回積分すると速度に、もう一回積分すると位置になります。
この積分は、上の数値積分の手法がそのまま使えます。
ただし、加速度が時々刻々と変化して速度や位置が変化する一方で、その時々刻々の位置、速度が力に影響する場合は、まとめてsumで計算する訳にはいかないので、積算していくことにします。
(たとえば、バネにおもりをつけて引っ張って離すと、バネののびに応じて力がかかり、その力で移動し、移動した位置に応じて力がかかる、というように、力と位置に相互の関係があります)
 |
| 抵抗のない自由落下 |
自由落下の式はおなじみの
です。これは、本来、
から時刻tで２回積分すると、
となって得られる式です。
これを試してみましょう。
まず、数式を書き換えます。
計算の間隔をΔtとして、Δtの間の速度でxを増やし、同様にvをaで、aは固定です。
| A | B | C | D | E | |
| 1 | Δt | 0.01 | |||
| 2 | 重力g | 9.8 | |||
| 3 | |||||
| 4 | i | 時刻 t[s] | 落下位置 x[m] | 落下速度 v[m/s] | 加速度 a[m/ss] |
| 5 | 0 | =A5*$C$1 | 0 | 0 | =$C$2 |
| 6 | =A5+1 | ↓ | =C5+D5*$C$1 | =D5+E5*$C$1 | ↓ |
| 7 | ↓ | ↓ | ↓ | ↓ | ↓(300行くらい) |
その上で、それぞれの値を求める式を作って、あとは「必要な数だけ」(必要な時間÷Δt、様子を見ながらxが○○になるまでなど)ならべます。
ここで、x0,v0に相当する、C5,D5はただ「０」を書いてあるわけではありません。これは、「シミュレーション開始時点で位置がゼロ、速度もゼロとする」というシミュレーションの開始条件です。
たとえば、C5に「１０」と書けば、位置が１０のところからスタートしますし、今回は座標軸が下向きなので、D5(v0)に「−５」と書けば上向きに５の速度でスタート(下に-5)するので、投げ上げになります。
数表だけではわかりにくいため、グラフにして見やすくします。
シミュレーションによってはこの段階も重要で「可視化」といいます。結果があっても、膨大な数字のかたまりではわかりづらいので。
たんに、B〜E列で散布図のグラフにしています。(点付きのグラフを選ぶと表示が大変なので線のみがよい)
 |
| 空気抵抗のある自由落下 |
物体の断面積に比例、速度の２乗にも比例
という性質があります。そのため、落下中の物体にはおおざっぱには、
(Rは定数、A断面)という式で表されます。Rは空気中とか、物体の形状(球とか板とか棒とか)で決まります。
鉄の玉が速く落ち、ピンポン球が遅く落ちる、というのは、鉄のほうが重く、(/m)だけ小さくなり、本来の重力加速度に近い速度で落ちやすいためです。 同じ重さでも、中空の玉とぎっしり詰まった玉では後者のほうが小さく、Aが小さくなり、速く落ちます。
真空中では空気抵抗なく、R=0になるので、等しく落ちるわけです。
この式は簡単には積分できません。積分した結果のvがaの式に入っています。
こういうときには数値計算が便利です。
| A | B | C | D | E | |
| 1 | Δt | 0.1 | |||
| 2 | 重力g | 9.8 | |||
| 3 | 抵抗 | 0.01 | |||
| 4 | i | 時刻 t[s] | 落下位置 x[m] | 落下速度 v[m/s] | 加速度 a[m/ss] |
| 5 | 0 | =A5*$C$1 | 0 | 0 | =$C$2-$C$3*D5*D5 |
| 6 | =A5+1 | ↓ | =C5+D5*$C$1 | =D5+E5*$C$1 | ↓ |
| 7 | ↓ | ↓ | ↓ | ↓ | ↓ |
※比較的時間がかかるため、Δtを0.05に変えています。誤差は大きくなりますが傾向は見えます。
結果を見ると、速度があるところで一定になります。これはmgと空気抵抗が釣り合ったところで加速しなくなるためです。
そうでないと、雨もとんでもない速度で降ってくることになります(笑)。
ここで、C3の抵抗(上式の(RA/m)相当)を変えると、グラフの形が変わります。０にすれば抵抗なしで先ほどの式と同じ、位置も放物線になります。
大きくすればするほど、早い段階で速度が一定になるようになります。
このように、ちょっとした数字の変更でいろいろな状況を作り出せるところもシミュレーションの魅力です。
ただ、「現実にはあり得ない状況」を作り出して満足してしまうこともあるのが怖いところで、常に気をつける必要があります。
バネとダンパ
 |
| バネとダンパによる振動 |
 |
| バネとダンパによる振動 |
 |
| バネによる振動(ダンパなし) |
 |
| 数式を少し改良 |
運動を決める方程式は、
です。バネは伸び(縮み)に比例した力で位置を戻そうとします。それにくわえてダンパ(車やバイクの車輪のところについている筒状、ピストン状のもの)は、速度に比例してブレーキをかけます。
| A | B | C | D | E | |
| 1 | Δt | 0.01 | |||
| 2 | m | 1 | |||
| 3 | k | 1 | |||
| 4 | c | 0.5 | |||
| 5 | i | 時刻t[s] | 位置 x[m] | 速度 v[m/s] | 加速度 a[m/ss] |
| 6 | 0 | =A6*$C$1 | 1 | 0 | =-($C$3/$C$2)*C6-($C$4/$C$2)*D6 |
| 7 | =A6+1 | ↓ | =C6+D6*$C$1 | =D6+E6*$C$1 | ↓ |
| 8 | ↓ | ↓ | ↓ | ↓ | ↓(1000行くらい) |
グラフをみると、揺れが徐々に収まっていきますが、これがダンパの効果です。
さて、ここで、c=0、ありがちなバネの振動にしてみましょう。
そのとき、振動の周期は
となることを高校あたりで習っているかもしれません。で、実際グラフを見るとそうなります。
しかし、妙なところがあります。青線で表されている「位置」が、最初の１より若干大きくなっています。
常識の範囲では、長く振動し続けることはあっても、揺れが大きくなるはずはありません。
これはまさにシミュレーションの誤差です。誤差の要因は計算の誤差もありますが、一般には使った数式の問題です。
今回の例では、加速度→速度→位置と順に計算されるため、加速度から位置まで２Δｔの時間差が生じ、悪影響が出ています。
そのため、Δtを小さくすれば誤差は減り、Δtを大きくすると誤差が大きくなります(試しに、Δt=0.1とかにすると大変なことに)。
ただ、Δtを小さくすることは計算の回数が増えることを意味し、特に表計算では限度があります。
もう一つの解決策として、数式の修正があります。たとえば、
| A | B | C | D | E | |
| 6 | 0 | =A6*$C$1 | 1 | 0 | =-($C$3/$C$2)*C6-($C$4/$C$2)*D6 |
| 7 | =A6+1 | ↓ | =C6+D7*$C$1 | =D6+E6*$C$1 | ↓ |
| 8 | ↓ | ↓ | ↓ | ↓ | ↓(1000行くらい) |
 |
| バネと摩擦による減衰振動 |
 |
| 摩擦のある場合 |
動摩擦は動く向きと逆方向に働きますので、
となります。ここで摩擦F=μmg (μ：摩擦係数)です。
この条件分けの部分を作るには、
IF(条件式,条件成立の時の数式,条件不成立の時の数式)
を使います。今回は、
=-(k/m)x+IF(v＞0,-μg,+μg)
に相当する数式を加速度に指定すればいいわけです。
ただ、こういう怪しげな条件を加えるほど、シミュレーションはどっかにいきやすくなります。Δt=0.01では、μを0.1程度にすると、得体の知れないグラフになります(右図はμ=0.01)。
今回の例では「静止摩擦」が入っていないことが一つの要因です。普通は、vの負号が変わるとき一度は停止するため、そこで静止摩擦に引っかかって最後は止まりますが、この式ではそれがありません。
入れてみるといいでしょう。IFの「(不)成立時の数式」に、またIFを使った数式もいれられます。
ただし、条件にv=0は使えません。計算誤差もありますが、ちょうどv=0になる時刻で計算するとは限りません。
そのため、たとえば、「|v|＜小さな数字」などの条件を使います。
Excelでは「abs(v)＜0.001」などにの条件式となります。
表計算ソフトを使って思考する
コンピュータの得意なことは計算です。
ですが、計算の仕方は教えなければなりません。
ここまでやってきたような計算は、かなり教えるべきことは少ない計算ですが、「○○の条件になる数値」、たとえば f(x)=0になるようなx(方程式の解)を求めようと思ったとき、その求め方を教えるのが難しいことがあります。
そういうときは、計算はコンピュータに任せ、直感的な判断をするのは人間がやる、という分業をする手があります。
基本的に、そういうばあい「こうだ」という明確な方法はあまりありません。
明確な方法があれば「教えられます」。
ここでは、実際の例をやってみます。
方程式の解を求めてみる
ある条件を満たす変数を探す、まさに方程式の解です。
なにか設計したりするときに、そういう計算が必要なことがあります。
ここでは、単純な３次方程式を題材に、解を求めてみます。
自分でやってみるための実習用方程式。例によって学生番号。
 |
| とりあえずグラフを書いてみる |
 |
| それらしい範囲を拡大する |
 |
| 細かく探していく |
方程式をf(x)=0とみます。この解は、グラフy=f(x)がy=0(x軸)と交わるところのxです。
で、まずはグラフを書いておおざっぱに見当をつけます。
見つけるべき特徴は、
ｙがプラスからマイナス、−から＋に変わるところ
です。その間にy=0があるはずだからです。
グラフはとりあえず、xを適当な範囲で増やしつつ、対応するyを計算するようにして、グラフ化します。
| A | B | |
| 1 | x | y |
| 2 | -10 | =A2*A2*A2-???*A2*A2+???*A2-??? |
| 3 | =A2+1 | ↓ |
もちろん、適当に増加する値(いままでiと表記)からxを作ってもＯＫですし、絶対参照で初期値や増分を自由に変えられるようにしておくと、使い道が広がります。
最終的に、プラスとマイナスが変わるところを見つけます。
右の例では、xが3.9〜4.2, 5.1〜5.4, 8.1〜8.4の３カ所です。
ステップ２：細かく探していく。
とりあえず、3.9と4.2の間にあることはわかっているので、１桁ずつ見つけていきましょう。
さっきのB列の数式をそのままコピーしてきて、A列の値を手動でいれて、Bが0に近づくようにします。
ここでは、マイナス→プラスなので、yがマイナスだったらその上、yがプラスだったらxを大きくしすぎ、ということになります。
（プラス→マイナスだと逆傾向）
1: 3.9〜4.2では4.1〜4.2らしい
2: 4.12ちょいくらい？
3: 4.123でyが1.3×10の-7乗とかなり小さいので、これでいいことにしてみる。
どこまで繰り返すかは、どこまで正確な数字が必要か、ということになります。
バカらしい方法に見えますが、とりあえず値が必要、というときはわりと便利な方法です。
二分法
 |
| 二分法でコンピュータに計算させる |
そこで、以下の方法を使います。
- 下限をxd,上限をxuとする。f(xd)×f(xu)＜0のはず(＋→−か−→＋なので、かける負になる)。
- 平均値xm=(xd+xu)/2を計算する。
- もし、f(xd)f(xm)＜0なら、xd〜xmを次に調べる。もし、f(xu)f(xm)＜0ならxm〜xuを調べる
- 以上を繰り返す
| A | B | C | D | E | F | |
| 1 | xd | y=f(xd) | xd | y=f(xd) | xd | y=f(xd) |
| 2 | xd初期 | =A2*A2*A2.. | =(A2+E2)/2 | =C2*C2*C2.. | xu初期 | =E2*E2*E2.. |
| 3 | =IF(B2*D2＜0,A2,C2) | ↓ | ↓ | ↓ | =IF(D2*F2＜0,E2,C2) | ↓ |
| 4 | ↓ | ↓ | ↓ | ↓ | ↓ | ↓ |
「=IF(B2*D2＜0,A2,C2)」は、もし、「B2*D2＜0ならば、A2の値を使う、そうでなければC2の値を使う」、という意味です。
「B2*D2＜0」はxdとxuの間をxm区切ってみて、そのxd側に答えがある、ということを意味します。そのため、次のチェック範囲はxdとxmの間なので、そのままxdを使います。一方、「そうでなければ」の条件は、チェック範囲がxmとxdの範囲になります。 つまり、いまのxdは捨てて、今度はxmをxdの代わりに使って、同じチェックをすればいいわけです。なので、xmであるC2を持ってきます。
この二分法は、１段の計算で、xの範囲は1/2に縮みます。10回もやれば、1/1000に、20回やれば１００万分の１までいきます。
また、最初のxu,xdを設定し直せば一瞬で他の区間の答えも出ます。
手でxを調整するよりは、じつはこれのほうが効率がいいはずですが、上の作表を思いつくのに時間がかかりそうなら、手動調整のほうが速いのも確かです。
ちなみに、区間内に複数の答えがあった場合には、適当に１個が見つかります。その意味で、ステップ１のグラフであたりをつける、が重要になります。
また、xd,xm,xuのうち２つは次段に引き継ぐので、一緒に対応するf(xd),f(xm),f(xu)も引き継ぐと、f(x)が複雑な時などにはスピードアップが期待できます。
時と場合で手段を切り替えることもまた、コンピュータを使いこなす秘訣です。
ロボットの挙動を確認する
 |
| 車輪移動ロボット |
ここでは、ロボットの車輪の回転速度を時々刻々と与えると、ロボットの挙動が出てくるように、計算式を並べています。
この計算をつかうと、本物のロボットの車輪の時々刻々のデータから、ロボットが今どこにいるかを推定することもできますし、ロボットを動かしてみなくても、大まかな動きを検討することができます。
ファイルを用意しておきますので、興味があれば、試してみるといいでしょう。
車輪ロボシミュレータ (2007/10/26, 66,048 bytes)
使い方：C列D列に左右の車軸の回転速度を入れます。
レポート課題
 |
| 振り子 |
 |
| 振り子シミュレーション例 |
振り子は、以下の解説に示すように、振り子の向いている角度θ(真下で0,真横で±π/2[rad])について、
または、
（単に角速度ωを用意して、分解しただけ）
という関係が成り立ちます。
この振り子の運動をシミュレートし、みてわかるような結果でグラフ化してください。
なお、振り子の初期角度は「80＋学籍番号下一桁[deg]」(80〜89)、振り子の長さは「3＋学籍番号百の位」(4or5)でつくってください。
※角度はdegからradに変換しないといけない。＜θは必ず[rad]
課題番号r03、ファイルはxls形式、"r03_学生番号.xls"にて、 期限は１１月９日 。
講義中に出た質問から
- d2θ/dt2って？ →θを２回tで微分するとこれになります。逆に、これを２回積分すると、θになります。
バネの例で、a→v→xとした過程と同じく、d2θ/dt2(角加速度)→dθ/dt(角速度,ω)→θ(角度)になります。
つまり、バネの例と積分する部分の計算はなにもかわりません。 - DEG→RADの仕方は？ ほぼ全員、１回目のレポートで、自分の手でやっているはずですが？(笑)。 斜方投射のあたり、見直してみましょう。/180*Pi()でできます。
- 振り子の長さや初期角度は簡単に変えられるようにする。
- θが小さいときは、sin(√(g/l)t)のただの正弦波で表すことができるとされている。(sinθをθと近似できる場合) それとの波形的な比較。
- 初期角度θ0と周期の関係を集めて(一種の簡単な実験)、グラフにしてみる。
- θ-tのグラフだけではなく、振り子の運動そのものに着目して、横方向xや上下動yなどのグラフを用意してみる。
- エネルギー保存則が成り立っているか、計算してみる。
- 空気抵抗が働いてもいいかも。
- そのほか、「私はここまでできます」と主張するようなネタ。
振り子の方程式
振り子が角度θ傾いているとします。
このとき、振り子には、θを減らす方向に力mgsinθがかかります。
振り子のその瞬間の速度をvとすると、f=maに当てはめて、
振り子は半径が振り子の長さlの円周状を動くため、その速度は傾く速度(角速度)ω×lとなります(円周の長さ＝半径×角度、を両辺とも時間で微分)。
ついでにmもとってしまうと
よって、
今回のレポートでは使いませんが、問題を簡単にする場合は、θが十分小さければsinθとθがほとんど同じとして、
と近似します。その場合、たとえば、
がこの式を満たします。
このページここで終了。
熊谷正朗 [→連絡]
東北学院大学 工学部 機械知能工学科 RDE ［｜ ］