#author("2016-10-26T08:03:08+09:00","","")
&size(24){表計算ソフトの応用その２};
このページについて
まえのページでの基本的操作をもとに、実際に工学部＆機械で役に立ちそうな計算方法、活用方法を紹介し、実際に試してみます。
#contents
**積算/差分/数値積分 [#s0cfe142]
ここでは、主に「積分」に関わるような計算を実例を通してみていきます。
***積算 [#ib299caf]
「道のり」＝「速さ」×「時間」
は、小学校あたりで習う、日常的な計算です。
実際のものの動きは、時々刻々と速度が変わります。
そこで、「速度○○で○○秒移動」という記録が多数あったときの計算を行ってみます。
速度の積算&br;
&ref(thumbFH_comp3_sumup1.gif);&br;
速度の積算&br;
&br;
速度を積算したイメージ図&br;
&ref(thumbFH_comp3_sumup2.gif);&br;
速度を積算したイメージ図&br;
&br;
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キロ弱)。
***差分 [#r2757c75]
差分&br;
&ref(thumbFH_comp3_sumup3.gif);&br;
差分&br;
&br;
逆に、位置がわかっている状態から速度変化を求めることもできます。
A B C D E
1 時刻 位置 移動量 経過時間 速度
2 時刻１ 位置１
3 時刻２ 位置２ =B3-B2 =A3-A2 =C3/D3
4 ： ： ↓ ↓ ↓
時刻とともに位置のデータが与えられているとき、単に、位置と時刻の差をとって、割るだけです。
これも、通常は「一定の時間間隔で計測」することが多く、その場合は経過時間が一定となるため、割り算も定数ですみます。
この計算もまた、計測制御の分野で実際によく使われる演算で、微分の代わりにも用いられます。
なお、差や速度の結果である、Ｃ〜Ｅ列で、２行目があけてあります。
これは趣味の問題で、３行目以降を上に詰めてもかまいません。 こういう計算をするときには、「時間が下に向かって進む」ため、「左斜め下＝気分的に時間が後ろ＝未来」の値を使うような印象があるため、このような形にしています。
***数値積分 [#y6efeeed]
先ほどの積算の計算をする際、積算すべき対象を自分で関数で用意することで、その関数の定積分(に近い値)を計算することができます。
この手法を数値積分、積分の数値(計算)解と呼びます。 一方、微分積分学の講義で習うような、数式のやりくりだけで、最後に値を入れれば結果が出るような結果を解析解といいます。
基本的に解析解は数式変形上の正しさを保って行うため、値も正確です。 しかし、数値解は普通は後述のように誤差が生じます。
誤差は出ますが、解析解のように悩まなくても、そもそも解析解がない場合でも、そこそこ答えが得られることが最大の利点です。
実際に、ここでは積分してみましょう。
y=2xの数値積分(0.1step)&br;
&ref(thumbFH_comp3_int10.gif);&br;
y=2xの数値積分(0.1step)&br;
y=2xの数値積分(0.1step)&br;
&ref(thumbFH_comp3_int1.gif);&br;
y=2xの数値積分(0.1step)&br;
y=2xの数値積分(0.01step)&br;
&ref(thumbFH_comp3_int2.gif);&br;
y=2xの数値積分(0.01step)&br;
まずは、簡単な例として、&br;
&ref(spreadsheet2_bhtml_eqn21.gif);&br;
を積分していってみます。
答えはいうまでもなく、&br;
&ref(spreadsheet2_bhtml_eqn22.gif);&br;
です(積分定数＝０)。
ここでは、先の例での時間にあたるｘを一定の幅で増加させ、そのときのｙを数式で与えます。 それ以外はいっしょです。
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までの計算値になる。
数値積分の誤差&br;
&ref(thumbFH_numint1.gif);&br;
数値積分の誤差&br;
この誤差の原因をもう少しよく考えてみます。
そもそも積分のおおざっぱな解釈は、「関数ｙ＝ｆ(ｘ)と、ｘ軸に囲まれた範囲の面積」です(y＞0ならそのもの)。
図は「解析値＝曲線の下側の面積＝ピンク」と「数値積分＝短冊の面積の合計＝水色」を図にしたものです。
この図では、曲線が増加しているときには、水色はピンクに足りていません。つまり、少なめに結果が出てしまいます。
逆に、減少時には水色がはみ出しているので、合計値は多めになります。
アップダウンの激しい関数の場合は、適当に帳消しになることが期待されますが、単調増加の場合には不足します。
ここで、短冊の横幅を細くしてみます。明らかにピンクの不足領域も、水色のはみ出しも小さくなっています。
これがΔｘを小さくしたときに誤差が減る理由です。
数値積分の誤差と計算方法&br;
&ref(thumbFH_numint2.gif);&br;
数値積分の誤差と計算方法&br;
この誤差は計算方法を変えると傾向が変わります。
図の一番上は「短冊の左端のｘで関数値を使う」です。この場合、先ほど同様、増加時に不足、減少時に過剰となります。
２番目の図は「短冊の右端のｘで関数値を使う」で、逆に、増加時に過剰、減少時に不足します。のでこの二つに大差はありません
３番目の図は「短冊の中点の値を使ってみた」です。なんとなくマシになる感じですが、頂点のところで誤差が出やすくなるなど、抜本的には改善しません。
４番目の図は「短冊をやめて台形にした」です。これは明らかに効果がありそうです。
実際にその効果が認められて「台形積分」と呼ばれる計算方法です。
計算量が増えそうですが、「１本目の短冊の右の高さ(y)」と「２本目の短冊の左の高さ」は同じ値です。ので、
積分＝（１左＋１右）Δｘ/2＋（２左＋２右）Δｘ/2＋..（ｎ左＋ｎ右）Δｘ/2
＝（Δｘ）｛１左/2＋１右(=２左)＋２右＋...(n-1)右＋ｎ右/2｝
となります。もとが、
積分＝（Δｘ）｛１左＋２左＋...＋ｎ左｝
なので、両端がすこし変わるだけです。途中の値がいらないなら、計算量がほとんど変わらないのに精度がよくなるという便利な方法です。
難点は、両端の例外部分(/2)を作り間違う可能性がある、というところです。
いまどきはコンピュータの速度が劇的に速いので、力任せにΔｘを細かくしてすませる、も一つの手です。
xsinx 0〜π の積分&br;
&ref(thumbFH_comp3_int3.gif);&br;
xsinx 0〜π の積分&br;
次にもう少し複雑な例として、&br;
&ref(spreadsheet2_bhtml_eqn23.gif);&br;
です。やれば、高校クラスの数学知識でも答えはでますが、ぱっとみて「積分したくねー」と思う程度の積分です。
まず、解析解を求めて置きます。&br;
&ref(spreadsheet2_bhtml_eqn24.gif);&br;
経過の理解は置いておいて、結果はπ=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 ↓ ↓ ↓ ↓
これで、分割１００個にあたる、No=99のところでみてみます。
積算結果は3.1413。解析解では3.1416なので、まあまああってます。
この分割を細かくしていけば、当然精度はよくなりますが、表計算の上で手頃にできるのは１００くらい、がんばっても実用になるのは１０００分割くらいです。
それ以上細かくして計算精度を上げるには、計算法を見直す(台形積分をつかってみるとか)ことも一つですが、プログラムを書けば、このくらいの積分なら100万分割もわけなく計算できます。これはプログラミングのところで今後やることになります。
この計算は先に述べたように、解析解を求めるのが面倒/不可能(自分には、を含む)なときにそれらしい値を出してくれます。
加えて、「解析解を出せたけど、正しいかが疑わしい」という場合に、検算として使えますし、予めこのあたりの値と目安をつけておくと解析しやすい場合もあります。
そういう「セカンドオピニオン」としての数値解、覚えておくと便利です。
すくなくとも、私の場合は、大学時代には数学の宿題や、研究で必要な数式の確認のため、よく、この手を使っていました。
確実に、着実に結果を求めていく必要がある場合、「常に２系統以上の答えの出し方を考え、結果の一致を確認する」ことは非常に重要です。
余裕があればぜひやってみましょう
•分割１０００にしてみる。
→$C$2を/1000書き換え、数式を１０００行くらいまで増やすだけ
•y=xcosxを積分してみる。
→C列の数式をさしかえ。解析解は２？&br;
&ref(spreadsheet2_bhtml_eqn25.gif);&br;
を−１〜１で積分してみる。
xを−１〜１になるように＆ｙ書き換え。√aはSqrt(a)で求まる(=Sqrt(1*1-B?*B?))。
結果はπ/2=1.57くらいになるはず。これ、何の積分？＜x,yでグラフを書けば一発でわかる？。
•台形積分をやってみる。
**シミュレーション [#f071872a]
***シミュレーションとは [#r67bdfc6]
機械におけるコンピュータの応用事例で多いものにコンピュータシミュレーションがあります。
シミュレーションとは、検討を加えたい対象の挙動を決める各種方程式を元に、順次計算を行っていき、その結果を評価するものです。 これにより、実物での実験をすることなく、様々な結果を得ることができるほか、実物では計測し得ないデータを得ることもできます。
たとえば、
•流体シミュレーション：自動車や飛行機の周りの空気の流れなど
•機械的シミュレーション：材料の変形、振動など。加工の模擬なども
•機構シミュレーション：メカの動きを計算、部品の衝突などが起きないかどうか
•ロボット制御シミュレーション：ロボットを模擬し、それを制御するにはどうするかを検討
•回路シミュレーション：設計した電子回路が適切に動作するかを検証
などがあります。いずれにしても、 •対象を的確に表す数式
•動作の条件(多くの場合、計算する領域の端をどうするかが重要で、境界条件と呼ばれる)
•計算の細かさ
などが重要項目で、これらに問題があると本来確認したいはずの現実と、シミュレーション結果がずれます。 時には、ちょっとのずれが全く異なる結果をもたらすこともあり、「シミュレーション結果によれば大丈夫」という場合には、「シミュレーションそのものが大丈夫」であるかの検証が必要です。
個人的な感想からすると、流体シミュレーションはかなり現実に近い結果が出るようですが、ロボット系、特に歩行ロボットは地面との断続的な接触部分をどうするか、滑りをどうするかが完璧には計算できず、往々にして実物と違う結果がでるようです。
ここでは、シミュレーションを実例で確認してみます。
***空気中の自由落下 [#xceb506f]
物体の運動は、&br;
&ref(spreadsheet2_bhtml_eqn101.gif);&br;
という運動方程式で表されます。言葉で書き直すと、
物体にかかる力＝質量×加速度
です。加速度は、速度の瞬間的な時間変化、つまり位置を時間で２回微分したものになります。
逆にみると、&br;
&ref(spreadsheet2_bhtml_eqn102.gif);&br;
加速度は力を質量で割ったもの、となります。
ものの運動のシミュレーションの基本はこの式になります。
時々刻々と物体にかかる力(重力とか駆動力とか)を求め(or入力し)、加速度にします。 これを１回積分すると速度に、もう一回積分すると位置になります。
この積分は、上の数値積分の手法がそのまま使えます。 ただし、加速度が時々刻々と変化して速度や位置が変化する一方で、その時々刻々の位置、速度が力に影響する場合は、まとめてsumで計算する訳にはいかないので、積算していくことにします。
(たとえば、バネにおもりをつけて引っ張って離すと、バネののびに応じて力がかかり、その力で移動し、移動した位置に応じて力がかかる、というように、力と位置に相互の関係があります)
抵抗のない自由落下&br;
&ref(thumbFH_comp4_11.gif);&br;
抵抗のない自由落下&br;
まず、簡単な自由落下をしてみます。
自由落下の式はおなじみの&br;
&ref(spreadsheet2.bhtml.eqn103.gif);&br;
です。これは、本来、&br;
&ref(spreadsheet2_bhtml_eqn104.gif);&br;
から時刻tで２回積分すると、&br;
&ref(spreadsheet2_bhtml_eqn105.gif);&br;
となって得られる式です。
これを試してみましょう。
まず、数式を書き換えます。&br;
&ref(spreadsheet2_bhtml_eqn106.gif);&br;
計算の間隔をΔ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列で散布図のグラフにしています。(点付きのグラフを選ぶと表示が大変なので線のみがよい)
空気抵抗のある自由落下&br;
&ref(comp4_12.gif);&br;
空気抵抗のある自由落下&br;
さて、現実には空気抵抗があります。空気抵抗は主に、
物体の断面積に比例、速度の２乗にも比例
という性質があります。そのため、落下中の物体にはおおざっぱには、&br;
&ref(spreadsheet2_bhtml_eqn107.gif);&br;
(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 ↓ ↓ ↓ ↓ ↓
赤で示したところが変更箇所です。E5の式を書き換えたら、忘れずにしたまでコピーしておきます。
※比較的時間がかかるため、Δtを0.05に変えています。誤差は大きくなりますが傾向は見えます。
結果を見ると、速度があるところで一定になります。これはmgと空気抵抗が釣り合ったところで加速しなくなるためです。 そうでないと、雨もとんでもない速度で降ってくることになります(笑)。 ここで、C3の抵抗(上式の(RA/m)相当)を変えると、グラフの形が変わります。０にすれば抵抗なしで先ほどの式と同じ、位置も放物線になります。 大きくすればするほど、早い段階で速度が一定になるようになります。
このように、ちょっとした数字の変更でいろいろな状況を作り出せるところもシミュレーションの魅力です。
ただ、「現実にはあり得ない状況」を作り出して満足してしまうこともあるのが怖いところで、常に気をつける必要があります。
**バネとダンパ [#aee839ff]
バネとダンパによる振動&br;
&ref(sim_sdm1.gif);&br;
バネとダンパによる振動&br;
バネとダンパによる振動&br;
&ref(thumbFH_comp4_20.gif);
バネとダンパによる振動%&br;
バネによる振動(ダンパなし)&br;
&ref(thumbFH_comp4_21.gif);
バネによる振動(ダンパなし)&br;
数式を少し改良&br;
&ref(thumbFH_comp4_22.gif);&br;
数式を少し改良&br;
二つ目の例として、バネとダンパという、機械ではおなじみのものをみてみましょう。
運動を決める方程式は、&br;
&ref(spreadsheet2_bhtml_eqn108.gif);&br;
です。バネは伸び(縮み)に比例した力で位置を戻そうとします。それにくわえてダンパ(車やバイクの車輪のところについている筒状、ピストン状のもの)は、速度に比例してブレーキをかけます。
&br;
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行くらい)
先ほどと変わったところは、基本的には、E列の加速度のみです。C6,D6はやはり初期位置、速度です(バネを1[m]引いて静止した状態でスタート)。
グラフをみると、揺れが徐々に収まっていきますが、これがダンパの効果です。
さて、ここで、c=0、ありがちなバネの振動にしてみましょう。
そのとき、振動の周期は
2\pi\sqrt{m/k}
となることを高校あたりで習っているかもしれません。で、実際グラフを見るとそうなります。
しかし、妙なところがあります。青線で表されている「位置」が、最初の１より若干大きくなっています。
常識の範囲では、長く振動し続けることはあっても、揺れが大きくなるはずはありません。
これはまさにシミュレーションの誤差です。誤差の要因は計算の誤差もありますが、一般には使った数式の問題です。 今回の例では、加速度→速度→位置と順に計算されるため、加速度から位置まで２Δｔの時間差が生じ、悪影響が出ています。
そのため、Δ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行くらい)
位置を求める式を直前の速度ではなく今の速度に代用します。 すると、先ほどの増加がなくなります(速度も今の加速度にするとx,v,aが互いに絡まる式になってExcelが混乱)。 見た目はＯＫですが、位置は厳密には「直前から現在までのΔｔのあいだの速度の積分」だけ増えるべきなので、「今の速度」にすればいい、というわけでもありません。ただ、誤差が少なくなって見える分だけ、無難とは言えます。
-----
バネと摩擦による減衰振動&br;
&ref(sim_sdm2.gif);&br;
バネと摩擦による減衰振動&br;
摩擦のある場合&br;
&ref(thumbFH_comp4_23.gif);&br;
摩擦のある場合&br;
発展的挑戦 ダンパの代わりに摩擦が働くことを考えます。
動摩擦は動く向きと逆方向に働きますので、&br;
&ref(spreadsheet2_bhtml_eqn110.gif);&br;
となります。ここで摩擦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」などにの条件式となります。
**表計算ソフトを使って思考する [#z6f987b0]
ここまではもっぱら指示だけして計算させまくってましたが、「コンピュータと一緒に」考えてみましょう。
コンピュータの得意なことは計算です。 ですが、計算の仕方は教えなければなりません。
ここまでやってきたような計算は、かなり教えるべきことは少ない計算ですが、「○○の条件になる数値」、たとえば f(x)=0になるようなx(方程式の解)を求めようと思ったとき、その求め方を教えるのが難しいことがあります。
そういうときは、計算はコンピュータに任せ、直感的な判断をするのは人間がやる、という分業をする手があります。
基本的に、そういうばあい「こうだ」という明確な方法はあまりありません。 明確な方法があれば「教えられます」。 ここでは、実際の例をやってみます。
**二分法 [#t8de93f7]
二分法でコンピュータに計算させる&br;
&ref(thumbFH_comp4_34.gif);&br;
二分法でコンピュータに計算させる&br;
先ほどのステップ２はある意味、行動パターンが決まっています。 「この間にある」とわかっているのxの範囲で、その間のxを試してみる、という形になります。 そのどちらかで、＋→＋や−→−だと、そこには解はなく、次に調べるのは＋→−か−→＋になる範囲です。
そこで、以下の方法を使います。
1.下限をxd,上限をxuとする。f(xd)×f(xu)＜0のはず(＋→−か−→＋なので、かける負になる)。
2.平均値xm=(xd+xu)/2を計算する。
3.もし、f(xd)f(xm)＜0なら、xd〜xmを次に調べる。もし、f(xu)f(xm)＜0ならxm〜xuを調べる
4.以上を繰り返す
この手順をExcelにすると、
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 ↓ ↓ ↓ ↓ ↓ ↓
となります。肝はA3およびE3のxd,xuの修正作業です。
「=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)が複雑な時などにはスピードアップが期待できます。
時と場合で手段を切り替えることもまた、コンピュータを使いこなす秘訣です。
**ロボットの挙動を確認する [#t17bb996]
車輪移動ロボット&br;
&ref(thumbFH_comp4_24.gif);&br;
車輪移動ロボット&br;
車輪が２個ついたようなロボットは、車輪が滑らないという前提ではその動きを計算で求めることができます。 →参考
ここでは、ロボットの車輪の回転速度を時々刻々と与えると、ロボットの挙動が出てくるように、計算式を並べています。
この計算をつかうと、本物のロボットの車輪の時々刻々のデータから、ロボットが今どこにいるかを推定することもできますし、ロボットを動かしてみなくても、大まかな動きを検討することができます。
ファイルを用意しておきますので、興味があれば、試してみるといいでしょう。
車輪ロボシミュレータ (2007/10/26, 66,048 bytes)
使い方：C列D列に左右の車軸の回転速度を入れます。
レポート課題
振り子&br;
&ref(thumbFH_sim_pend.gif);&br;
振り子 &br;
振り子シミュレーション例&br;
&ref(thumbFH_comp4_41.gif);&br;
振り子シミュレーション例&br;
表計算ソフトの部分のまとめとして、振り子のシミュレーションを課題とします。
振り子は、以下の解説に示すように、振り子の向いている角度θ(真下で0,真横で±π/2[rad])について、&br;
&ref(spreadsheet2_bhtml_eqn51.gif);&br;
または、&br;
&ref(spreadsheet2_bhtml_eqn52.gif);&br;
（単に角速度ωを用意して、分解しただけ）
という関係が成り立ちます。
この振り子の運動をシミュレートし、みてわかるような結果でグラフ化してください。
なお、振り子の初期角度は「80＋学籍番号下一桁[deg]」(80〜89)、振り子の長さは「3＋学籍番号百の位」(4or5)でつくってください。
※角度はdegからradに変換しないといけない。＜θは必ず[rad]
***質問からの考察 [#q3e7a4db]
•d2θ/dt2って？ →θを２回tで微分するとこれになります。逆に、これを２回積分すると、θになります。 バネの例で、a→v→xとした過程と同じく、d2θ/dt2(角加速度)→dθ/dt(角速度,ω)→θ(角度)になります。
つまり、バネの例と積分する部分の計算はなにもかわりません。
•DEG→RADの仕方は？ ほぼ全員、１回目のレポートで、自分の手でやっているはずですが？(笑)。 斜方投射のあたり、見直してみましょう。/180*Pi()でできます。
たとえば以下に示すような、各自の工夫を期待します(必須ではありませんが)。
•振り子の長さや初期角度は簡単に変えられるようにする。
•θが小さいときは、sin(√(g/l)t)のただの正弦波で表すことができるとされている。(sinθをθと近似できる場合) それとの波形的な比較。
•初期角度θ0と周期の関係を集めて(一種の簡単な実験)、グラフにしてみる。
•θ-tのグラフだけではなく、振り子の運動そのものに着目して、横方向xや上下動yなどのグラフを用意してみる。
•エネルギー保存則が成り立っているか、計算してみる。
•空気抵抗が働いてもいいかも。
このページここで終了。
&br;
//http://www.mech.tohoku-gakuin.ac.jp/rde/contents/course/computer/spreadsheet2.html
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
