表計算ソフト/表計算ソフトの応用その2
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
&size(24){表計算ソフトの応用その2};
このページについて
まえのページでの基本的操作をもとに、実際に工学部&機械で...
#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 速度1 時間1 =B2*C2 =D2 =C2
3 =A2+1 速度2 時間2 ↓ =E2+D3 =F2+C3
4 ↓ : : ↓ ↓ ↓
B列C列には「速度○○で○○秒移動」の実際のデータを順番に入...
D列は、その「○○秒間の移動距離」です(単に速度*時間)。
これをE列で累計していきます。もちろん、最終的なトータル...
表では「○○秒」のところがすべて2秒です。 この作表ならもち...
最初から「○○秒おき」がわかっていれば(それが普通)、わざわ...
ちなみに、この計算は当たり前なようで、計測や制御の分野な...
数学的には積分に近いため、後述のような数値演算に用いられ...
さて、実際にこの表のように、速度がすぱっと変わるかという...
おそらくは、右図の破線のような感じで速度が徐々に増え、一...
しかし、データはこれしかないので、とりあえず、水色の階段...
どうみてもずれがありますが、細かいことは後述します。
蛇足:
工学部だったら「秒速○○メートル」といったときに、「あ、あ...
ちなみに、歩く速度が 1[m/s]強(時速4キロ)です。世界最速の...
***差分 [#r2757c75]
差分&br;
&ref(thumbFH_comp3_sumup3.gif);&br;
差分&br;
&br;
逆に、位置がわかっている状態から速度変化を求めることもで...
A B C D E
1 時刻 位置 移動量 経過時間 速度
2 時刻1 位置1
3 時刻2 位置2 =B3-B2 =A3-A2 =C3/D3
4 : : ↓ ↓ ↓
時刻とともに位置のデータが与えられているとき、単に、位置...
これも、通常は「一定の時間間隔で計測」することが多く、そ...
この計算もまた、計測制御の分野で実際によく使われる演算で...
なお、差や速度の結果である、C〜E列で、2行目があけてあ...
これは趣味の問題で、3行目以降を上に詰めてもかまいません...
***数値積分 [#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;
です(積分定数=0)。
ここでは、先の例での時間にあたるxを一定の幅で増加させ、...
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 ↓ ↓ ↓ ↓ ↓ ↓
この作表でのポイントは以下の通りです。
•xの増加分Δxを特定のセルに置き、絶対参照($C$2)で使うよ...
→Δxをあれこれ変えて結果を確認できる。
•xもΔxをもとに作るようにする。B5=0として下向きにΔxを加...
•積分したいyをx(B列)の数式で作る。
•yΔxをつくる。
•E列で積算している。計算上はE6=E5+C6*$C$2のように、D列...
•補足:F列は比較のために、本来あるべき(解析値)y=x*x...
•補足:グラフは積算値E列と解析値F列を表示。
•補足:1行目、3行目の空白はなんとなく見やすくするための...
比較すると、右2番目の図のグラフのように、微妙にずれがあ...
このずれが、数値積分をする際に避けられない誤差です。
単純には、Δxを小さくすれば、誤差は減ります(右図3番目)。
Δx=0.1のときはx=3.0までの積分に当たるx=2.9(※)のところ...
単純な傾向としては、Δxを10分の1に細かくすると誤差は概...
※x=0.0からx=3.0までを30等分した状態なので、30個目...
数値積分の誤差&br;
&ref(thumbFH_numint1.gif);&br;
数値積分の誤差&br;
この誤差の原因をもう少しよく考えてみます。
そもそも積分のおおざっぱな解釈は、「関数y=f(x)と、x...
図は「解析値=曲線の下側の面積=ピンク」と「数値積分=短...
この図では、曲線が増加しているときには、水色はピンクに足...
逆に、減少時には水色がはみ出しているので、合計値は多めに...
アップダウンの激しい関数の場合は、適当に帳消しになること...
ここで、短冊の横幅を細くしてみます。明らかにピンクの不足...
これがΔxを小さくしたときに誤差が減る理由です。
数値積分の誤差と計算方法&br;
&ref(thumbFH_numint2.gif);&br;
数値積分の誤差と計算方法&br;
この誤差は計算方法を変えると傾向が変わります。
図の一番上は「短冊の左端のxで関数値を使う」です。この場...
2番目の図は「短冊の右端のxで関数値を使う」で、逆に、増...
3番目の図は「短冊の中点の値を使ってみた」です。なんとな...
4番目の図は「短冊をやめて台形にした」です。これは明らか...
実際にその効果が認められて「台形積分」と呼ばれる計算方法...
計算量が増えそうですが、「1本目の短冊の右の高さ(y)」と「...
積分=(1左+1右)Δx/2+(2左+2右)Δx/2+..(n左...
=(Δx){1左/2+1右(=2左)+2右+...(n-1)右+n右/2}
となります。もとが、
積分=(Δx){1左+2左+...+n左}
なので、両端がすこし変わるだけです。途中の値がいらないな...
難点は、両端の例外部分(/2)を作り間違う可能性がある、とい...
いまどきはコンピュータの速度が劇的に速いので、力任せにΔx...
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 ↓ ↓ ↓ ↓
これで、分割100個にあたる、No=99のところでみてみます。
積算結果は3.1413。解析解では3.1416なので、まあまああって...
この分割を細かくしていけば、当然精度はよくなりますが、表...
それ以上細かくして計算精度を上げるには、計算法を見直す(台...
この計算は先に述べたように、解析解を求めるのが面倒/不可能...
加えて、「解析解を出せたけど、正しいかが疑わしい」という...
そういう「セカンドオピニオン」としての数値解、覚えておく...
すくなくとも、私の場合は、大学時代には数学の宿題や、研究...
確実に、着実に結果を求めていく必要がある場合、「常に2系...
余裕があればぜひやってみましょう
•分割1000にしてみる。
→$C$2を/1000書き換え、数式を1000行くらいまで増やすだけ
•y=xcosxを積分してみる。
→C列の数式をさしかえ。解析解は2?&br;
&ref(spreadsheet2_bhtml_eqn25.gif);&br;
を−1〜1で積分してみる。
xを−1〜1になるように&y書き換え。√aはSqrt(a)で求まる(=...
結果はπ/2=1.57くらいになるはず。これ、何の積分?<x,yでグ...
•台形積分をやってみる。
**シミュレーション [#f071872a]
***シミュレーションとは [#r67bdfc6]
機械におけるコンピュータの応用事例で多いものにコンピュー...
シミュレーションとは、検討を加えたい対象の挙動を決める各...
たとえば、
•流体シミュレーション:自動車や飛行機の周りの空気の流れな...
•機械的シミュレーション:材料の変形、振動など。加工の模擬...
•機構シミュレーション:メカの動きを計算、部品の衝突などが...
•ロボット制御シミュレーション:ロボットを模擬し、それを制...
•回路シミュレーション:設計した電子回路が適切に動作するか...
などがあります。いずれにしても、 •対象を的確に表す数式
•動作の条件(多くの場合、計算する領域の端をどうするかが重...
•計算の細かさ
などが重要項目で、これらに問題があると本来確認したいはず...
個人的な感想からすると、流体シミュレーションはかなり現実...
ここでは、シミュレーションを実例で確認してみます。
***空気中の自由落下 [#xceb506f]
物体の運動は、&br;
&ref(spreadsheet2_bhtml_eqn101.gif);&br;
という運動方程式で表されます。言葉で書き直すと、
物体にかかる力=質量×加速度
です。加速度は、速度の瞬間的な時間変化、つまり位置を時間...
逆にみると、&br;
&ref(spreadsheet2_bhtml_eqn102.gif);&br;
加速度は力を質量で割ったもの、となります。
ものの運動のシミュレーションの基本はこの式になります。
時々刻々と物体にかかる力(重力とか駆動力とか)を求め(or入力...
この積分は、上の数値積分の手法がそのまま使えます。 ただし...
(たとえば、バネにおもりをつけて引っ張って離すと、バネのの...
抵抗のない自由落下&br;
&ref(thumbFH_comp4_11.gif);&br;
抵抗のない自由落下&br;
まず、簡単な自由落下をしてみます。
自由落下の式はおなじみの&br;
&ref(spreadsheet2.bhtml.eqn103.gif);&br;
です。これは、本来、&br;
&ref(spreadsheet2_bhtml_eqn104.gif);&br;
から時刻tで2回積分すると、&br;
&ref(spreadsheet2_bhtml_eqn105.gif);&br;
となって得られる式です。
これを試してみましょう。
まず、数式を書き換えます。&br;
&ref(spreadsheet2_bhtml_eqn106.gif);&br;
計算の間隔をΔtとして、Δtの間の速度でxを増やし、同様にvを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行くらい)
表を作るときは、時間とともに変化する値を横に変数として並...
その上で、それぞれの値を求める式を作って、あとは「必要な...
ここで、x0,v0に相当する、C5,D5はただ「0」を書いてあるわ...
たとえば、C5に「10」と書けば、位置が10のところからス...
数表だけではわかりにくいため、グラフにして見やすくします。
シミュレーションによってはこの段階も重要で「可視化」とい...
たんに、B〜E列で散布図のグラフにしています。(点付きのグラ...
空気抵抗のある自由落下&br;
&ref(comp4_12.gif);&br;
空気抵抗のある自由落下&br;
さて、現実には空気抵抗があります。空気抵抗は主に、
物体の断面積に比例、速度の2乗にも比例
という性質があります。そのため、落下中の物体にはおおざっ...
&ref(spreadsheet2_bhtml_eqn107.gif);&br;
(Rは定数、A断面)
という式で表されます。Rは空気中とか、物体の形状(球とか板...
鉄の玉が速く落ち、ピンポン球が遅く落ちる、というのは、鉄...
真空中では空気抵抗なく、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...
このように、ちょっとした数字の変更でいろいろな状況を作り...
ただ、「現実にはあり得ない状況」を作り出して満足してしま...
**バネとダンパ [#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列の加速度のみです...
グラフをみると、揺れが徐々に収まっていきますが、これがダ...
さて、ここで、c=0、ありがちなバネの振動にしてみましょう。
そのとき、振動の周期は
2\pi\sqrt{m/k}
となることを高校あたりで習っているかもしれません。で、実...
しかし、妙なところがあります。青線で表されている「位置」...
常識の範囲では、長く振動し続けることはあっても、揺れが大...
これはまさにシミュレーションの誤差です。誤差の要因は計算...
そのため、Δtを小さくすれば誤差は減り、Δtを大きくすると誤...
ただ、Δ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行くらい)
位置を求める式を直前の速度ではなく今の速度に代用します。 ...
-----
バネと摩擦による減衰振動&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)
に相当する数式を加速度に指定すればいいわけです。
ただ、こういう怪しげな条件を加えるほど、シミュレーション...
今回の例では「静止摩擦」が入っていないことが一つの要因で...
入れてみるといいでしょう。IFの「(不)成立時の数式」に、ま...
ただし、条件にv=0は使えません。計算誤差もありますが、ちょ...
そのため、たとえば、「|v|<小さな数字」などの条件を使いま...
Excelでは「abs(v)<0.001」などにの条件式となります。
**表計算ソフトを使って思考する [#z6f987b0]
ここまではもっぱら指示だけして計算させまくってましたが、...
コンピュータの得意なことは計算です。 ですが、計算の仕方は...
ここまでやってきたような計算は、かなり教えるべきことは少...
そういうときは、計算はコンピュータに任せ、直感的な判断を...
基本的に、そういうばあい「こうだ」という明確な方法はあま...
**二分法 [#t8de93f7]
二分法でコンピュータに計算させる&br;
&ref(thumbFH_comp4_34.gif);&br;
二分法でコンピュータに計算させる&br;
先ほどのステップ2はある意味、行動パターンが決まっていま...
そこで、以下の方法を使います。
1.下限をxd,上限をxuとする。f(xd)×f(xu)<0のはず(+→−か−→...
2.平均値xm=(xd+xu)/2を計算する。
3.もし、f(xd)f(xm)<0なら、xd〜xmを次に調べる。もし、f(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...
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の...
「B2*D2<0」はxdとxuの間をxm区切ってみて、そのxd側に答え...
この二分法は、1段の計算で、xの範囲は1/2に縮みます。10回...
また、最初のxu,xdを設定し直せば一瞬で他の区間の答えも出ま...
手でxを調整するよりは、じつはこれのほうが効率がいいはず...
ちなみに、区間内に複数の答えがあった場合には、適当に1個...
また、xd,xm,xuのうち2つは次段に引き継ぐので、一緒に対応...
時と場合で手段を切り替えることもまた、コンピュータを使い...
**ロボットの挙動を確認する [#t17bb996]
車輪移動ロボット&br;
&ref(thumbFH_comp4_24.gif);&br;
車輪移動ロボット&br;
車輪が2個ついたようなロボットは、車輪が滑らないという前...
ここでは、ロボットの車輪の回転速度を時々刻々と与えると、...
この計算をつかうと、本物のロボットの車輪の時々刻々のデー...
ファイルを用意しておきますので、興味があれば、試してみる...
車輪ロボシミュレータ (2007/10/26, 66,048 bytes)
使い方:C列D列に左右の車軸の回転速度を入れます。
レポート課題
振り子&br;
&ref(thumbFH_sim_pend.gif);&br;
振り子 &br;
振り子シミュレーション例&br;
&ref(thumbFH_comp4_41.gif);&br;
振り子シミュレーション例&br;
表計算ソフトの部分のまとめとして、振り子のシミュレーショ...
振り子は、以下の解説に示すように、振り子の向いている角度θ...
&ref(spreadsheet2_bhtml_eqn51.gif);&br;
または、&br;
&ref(spreadsheet2_bhtml_eqn52.gif);&br;
(単に角速度ωを用意して、分解しただけ)
という関係が成り立ちます。
この振り子の運動をシミュレートし、みてわかるような結果で...
なお、振り子の初期角度は「80+学籍番号下一桁[deg]」(80〜8...
※角度はdegからradに変換しないといけない。<θは必ず[rad]
***質問からの考察 [#q3e7a4db]
•d2θ/dt2って? →θを2回tで微分するとこれになります。逆に...
つまり、バネの例と積分する部分の計算はなにもかわりません。
•DEG→RADの仕方は? ほぼ全員、1回目のレポートで、自分の...
たとえば以下に示すような、各自の工夫を期待します(必須では...
•振り子の長さや初期角度は簡単に変えられるようにする。
•θが小さいときは、sin(√(g/l)t)のただの正弦波で表すことが...
•初期角度θ0と周期の関係を集めて(一種の簡単な実験)、グラフ...
•θ-tのグラフだけではなく、振り子の運動そのものに着目して...
•エネルギー保存則が成り立っているか、計算してみる。
•空気抵抗が働いてもいいかも。
このページここで終了。
&br;
//http://www.mech.tohoku-gakuin.ac.jp/rde/contents/course...
終了行:
&size(24){表計算ソフトの応用その2};
このページについて
まえのページでの基本的操作をもとに、実際に工学部&機械で...
#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 速度1 時間1 =B2*C2 =D2 =C2
3 =A2+1 速度2 時間2 ↓ =E2+D3 =F2+C3
4 ↓ : : ↓ ↓ ↓
B列C列には「速度○○で○○秒移動」の実際のデータを順番に入...
D列は、その「○○秒間の移動距離」です(単に速度*時間)。
これをE列で累計していきます。もちろん、最終的なトータル...
表では「○○秒」のところがすべて2秒です。 この作表ならもち...
最初から「○○秒おき」がわかっていれば(それが普通)、わざわ...
ちなみに、この計算は当たり前なようで、計測や制御の分野な...
数学的には積分に近いため、後述のような数値演算に用いられ...
さて、実際にこの表のように、速度がすぱっと変わるかという...
おそらくは、右図の破線のような感じで速度が徐々に増え、一...
しかし、データはこれしかないので、とりあえず、水色の階段...
どうみてもずれがありますが、細かいことは後述します。
蛇足:
工学部だったら「秒速○○メートル」といったときに、「あ、あ...
ちなみに、歩く速度が 1[m/s]強(時速4キロ)です。世界最速の...
***差分 [#r2757c75]
差分&br;
&ref(thumbFH_comp3_sumup3.gif);&br;
差分&br;
&br;
逆に、位置がわかっている状態から速度変化を求めることもで...
A B C D E
1 時刻 位置 移動量 経過時間 速度
2 時刻1 位置1
3 時刻2 位置2 =B3-B2 =A3-A2 =C3/D3
4 : : ↓ ↓ ↓
時刻とともに位置のデータが与えられているとき、単に、位置...
これも、通常は「一定の時間間隔で計測」することが多く、そ...
この計算もまた、計測制御の分野で実際によく使われる演算で...
なお、差や速度の結果である、C〜E列で、2行目があけてあ...
これは趣味の問題で、3行目以降を上に詰めてもかまいません...
***数値積分 [#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;
です(積分定数=0)。
ここでは、先の例での時間にあたるxを一定の幅で増加させ、...
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 ↓ ↓ ↓ ↓ ↓ ↓
この作表でのポイントは以下の通りです。
•xの増加分Δxを特定のセルに置き、絶対参照($C$2)で使うよ...
→Δxをあれこれ変えて結果を確認できる。
•xもΔxをもとに作るようにする。B5=0として下向きにΔxを加...
•積分したいyをx(B列)の数式で作る。
•yΔxをつくる。
•E列で積算している。計算上はE6=E5+C6*$C$2のように、D列...
•補足:F列は比較のために、本来あるべき(解析値)y=x*x...
•補足:グラフは積算値E列と解析値F列を表示。
•補足:1行目、3行目の空白はなんとなく見やすくするための...
比較すると、右2番目の図のグラフのように、微妙にずれがあ...
このずれが、数値積分をする際に避けられない誤差です。
単純には、Δxを小さくすれば、誤差は減ります(右図3番目)。
Δx=0.1のときはx=3.0までの積分に当たるx=2.9(※)のところ...
単純な傾向としては、Δxを10分の1に細かくすると誤差は概...
※x=0.0からx=3.0までを30等分した状態なので、30個目...
数値積分の誤差&br;
&ref(thumbFH_numint1.gif);&br;
数値積分の誤差&br;
この誤差の原因をもう少しよく考えてみます。
そもそも積分のおおざっぱな解釈は、「関数y=f(x)と、x...
図は「解析値=曲線の下側の面積=ピンク」と「数値積分=短...
この図では、曲線が増加しているときには、水色はピンクに足...
逆に、減少時には水色がはみ出しているので、合計値は多めに...
アップダウンの激しい関数の場合は、適当に帳消しになること...
ここで、短冊の横幅を細くしてみます。明らかにピンクの不足...
これがΔxを小さくしたときに誤差が減る理由です。
数値積分の誤差と計算方法&br;
&ref(thumbFH_numint2.gif);&br;
数値積分の誤差と計算方法&br;
この誤差は計算方法を変えると傾向が変わります。
図の一番上は「短冊の左端のxで関数値を使う」です。この場...
2番目の図は「短冊の右端のxで関数値を使う」で、逆に、増...
3番目の図は「短冊の中点の値を使ってみた」です。なんとな...
4番目の図は「短冊をやめて台形にした」です。これは明らか...
実際にその効果が認められて「台形積分」と呼ばれる計算方法...
計算量が増えそうですが、「1本目の短冊の右の高さ(y)」と「...
積分=(1左+1右)Δx/2+(2左+2右)Δx/2+..(n左...
=(Δx){1左/2+1右(=2左)+2右+...(n-1)右+n右/2}
となります。もとが、
積分=(Δx){1左+2左+...+n左}
なので、両端がすこし変わるだけです。途中の値がいらないな...
難点は、両端の例外部分(/2)を作り間違う可能性がある、とい...
いまどきはコンピュータの速度が劇的に速いので、力任せにΔx...
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 ↓ ↓ ↓ ↓
これで、分割100個にあたる、No=99のところでみてみます。
積算結果は3.1413。解析解では3.1416なので、まあまああって...
この分割を細かくしていけば、当然精度はよくなりますが、表...
それ以上細かくして計算精度を上げるには、計算法を見直す(台...
この計算は先に述べたように、解析解を求めるのが面倒/不可能...
加えて、「解析解を出せたけど、正しいかが疑わしい」という...
そういう「セカンドオピニオン」としての数値解、覚えておく...
すくなくとも、私の場合は、大学時代には数学の宿題や、研究...
確実に、着実に結果を求めていく必要がある場合、「常に2系...
余裕があればぜひやってみましょう
•分割1000にしてみる。
→$C$2を/1000書き換え、数式を1000行くらいまで増やすだけ
•y=xcosxを積分してみる。
→C列の数式をさしかえ。解析解は2?&br;
&ref(spreadsheet2_bhtml_eqn25.gif);&br;
を−1〜1で積分してみる。
xを−1〜1になるように&y書き換え。√aはSqrt(a)で求まる(=...
結果はπ/2=1.57くらいになるはず。これ、何の積分?<x,yでグ...
•台形積分をやってみる。
**シミュレーション [#f071872a]
***シミュレーションとは [#r67bdfc6]
機械におけるコンピュータの応用事例で多いものにコンピュー...
シミュレーションとは、検討を加えたい対象の挙動を決める各...
たとえば、
•流体シミュレーション:自動車や飛行機の周りの空気の流れな...
•機械的シミュレーション:材料の変形、振動など。加工の模擬...
•機構シミュレーション:メカの動きを計算、部品の衝突などが...
•ロボット制御シミュレーション:ロボットを模擬し、それを制...
•回路シミュレーション:設計した電子回路が適切に動作するか...
などがあります。いずれにしても、 •対象を的確に表す数式
•動作の条件(多くの場合、計算する領域の端をどうするかが重...
•計算の細かさ
などが重要項目で、これらに問題があると本来確認したいはず...
個人的な感想からすると、流体シミュレーションはかなり現実...
ここでは、シミュレーションを実例で確認してみます。
***空気中の自由落下 [#xceb506f]
物体の運動は、&br;
&ref(spreadsheet2_bhtml_eqn101.gif);&br;
という運動方程式で表されます。言葉で書き直すと、
物体にかかる力=質量×加速度
です。加速度は、速度の瞬間的な時間変化、つまり位置を時間...
逆にみると、&br;
&ref(spreadsheet2_bhtml_eqn102.gif);&br;
加速度は力を質量で割ったもの、となります。
ものの運動のシミュレーションの基本はこの式になります。
時々刻々と物体にかかる力(重力とか駆動力とか)を求め(or入力...
この積分は、上の数値積分の手法がそのまま使えます。 ただし...
(たとえば、バネにおもりをつけて引っ張って離すと、バネのの...
抵抗のない自由落下&br;
&ref(thumbFH_comp4_11.gif);&br;
抵抗のない自由落下&br;
まず、簡単な自由落下をしてみます。
自由落下の式はおなじみの&br;
&ref(spreadsheet2.bhtml.eqn103.gif);&br;
です。これは、本来、&br;
&ref(spreadsheet2_bhtml_eqn104.gif);&br;
から時刻tで2回積分すると、&br;
&ref(spreadsheet2_bhtml_eqn105.gif);&br;
となって得られる式です。
これを試してみましょう。
まず、数式を書き換えます。&br;
&ref(spreadsheet2_bhtml_eqn106.gif);&br;
計算の間隔をΔtとして、Δtの間の速度でxを増やし、同様にvを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行くらい)
表を作るときは、時間とともに変化する値を横に変数として並...
その上で、それぞれの値を求める式を作って、あとは「必要な...
ここで、x0,v0に相当する、C5,D5はただ「0」を書いてあるわ...
たとえば、C5に「10」と書けば、位置が10のところからス...
数表だけではわかりにくいため、グラフにして見やすくします。
シミュレーションによってはこの段階も重要で「可視化」とい...
たんに、B〜E列で散布図のグラフにしています。(点付きのグラ...
空気抵抗のある自由落下&br;
&ref(comp4_12.gif);&br;
空気抵抗のある自由落下&br;
さて、現実には空気抵抗があります。空気抵抗は主に、
物体の断面積に比例、速度の2乗にも比例
という性質があります。そのため、落下中の物体にはおおざっ...
&ref(spreadsheet2_bhtml_eqn107.gif);&br;
(Rは定数、A断面)
という式で表されます。Rは空気中とか、物体の形状(球とか板...
鉄の玉が速く落ち、ピンポン球が遅く落ちる、というのは、鉄...
真空中では空気抵抗なく、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...
このように、ちょっとした数字の変更でいろいろな状況を作り...
ただ、「現実にはあり得ない状況」を作り出して満足してしま...
**バネとダンパ [#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列の加速度のみです...
グラフをみると、揺れが徐々に収まっていきますが、これがダ...
さて、ここで、c=0、ありがちなバネの振動にしてみましょう。
そのとき、振動の周期は
2\pi\sqrt{m/k}
となることを高校あたりで習っているかもしれません。で、実...
しかし、妙なところがあります。青線で表されている「位置」...
常識の範囲では、長く振動し続けることはあっても、揺れが大...
これはまさにシミュレーションの誤差です。誤差の要因は計算...
そのため、Δtを小さくすれば誤差は減り、Δtを大きくすると誤...
ただ、Δ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行くらい)
位置を求める式を直前の速度ではなく今の速度に代用します。 ...
-----
バネと摩擦による減衰振動&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)
に相当する数式を加速度に指定すればいいわけです。
ただ、こういう怪しげな条件を加えるほど、シミュレーション...
今回の例では「静止摩擦」が入っていないことが一つの要因で...
入れてみるといいでしょう。IFの「(不)成立時の数式」に、ま...
ただし、条件にv=0は使えません。計算誤差もありますが、ちょ...
そのため、たとえば、「|v|<小さな数字」などの条件を使いま...
Excelでは「abs(v)<0.001」などにの条件式となります。
**表計算ソフトを使って思考する [#z6f987b0]
ここまではもっぱら指示だけして計算させまくってましたが、...
コンピュータの得意なことは計算です。 ですが、計算の仕方は...
ここまでやってきたような計算は、かなり教えるべきことは少...
そういうときは、計算はコンピュータに任せ、直感的な判断を...
基本的に、そういうばあい「こうだ」という明確な方法はあま...
**二分法 [#t8de93f7]
二分法でコンピュータに計算させる&br;
&ref(thumbFH_comp4_34.gif);&br;
二分法でコンピュータに計算させる&br;
先ほどのステップ2はある意味、行動パターンが決まっていま...
そこで、以下の方法を使います。
1.下限をxd,上限をxuとする。f(xd)×f(xu)<0のはず(+→−か−→...
2.平均値xm=(xd+xu)/2を計算する。
3.もし、f(xd)f(xm)<0なら、xd〜xmを次に調べる。もし、f(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...
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の...
「B2*D2<0」はxdとxuの間をxm区切ってみて、そのxd側に答え...
この二分法は、1段の計算で、xの範囲は1/2に縮みます。10回...
また、最初のxu,xdを設定し直せば一瞬で他の区間の答えも出ま...
手でxを調整するよりは、じつはこれのほうが効率がいいはず...
ちなみに、区間内に複数の答えがあった場合には、適当に1個...
また、xd,xm,xuのうち2つは次段に引き継ぐので、一緒に対応...
時と場合で手段を切り替えることもまた、コンピュータを使い...
**ロボットの挙動を確認する [#t17bb996]
車輪移動ロボット&br;
&ref(thumbFH_comp4_24.gif);&br;
車輪移動ロボット&br;
車輪が2個ついたようなロボットは、車輪が滑らないという前...
ここでは、ロボットの車輪の回転速度を時々刻々と与えると、...
この計算をつかうと、本物のロボットの車輪の時々刻々のデー...
ファイルを用意しておきますので、興味があれば、試してみる...
車輪ロボシミュレータ (2007/10/26, 66,048 bytes)
使い方:C列D列に左右の車軸の回転速度を入れます。
レポート課題
振り子&br;
&ref(thumbFH_sim_pend.gif);&br;
振り子 &br;
振り子シミュレーション例&br;
&ref(thumbFH_comp4_41.gif);&br;
振り子シミュレーション例&br;
表計算ソフトの部分のまとめとして、振り子のシミュレーショ...
振り子は、以下の解説に示すように、振り子の向いている角度θ...
&ref(spreadsheet2_bhtml_eqn51.gif);&br;
または、&br;
&ref(spreadsheet2_bhtml_eqn52.gif);&br;
(単に角速度ωを用意して、分解しただけ)
という関係が成り立ちます。
この振り子の運動をシミュレートし、みてわかるような結果で...
なお、振り子の初期角度は「80+学籍番号下一桁[deg]」(80〜8...
※角度はdegからradに変換しないといけない。<θは必ず[rad]
***質問からの考察 [#q3e7a4db]
•d2θ/dt2って? →θを2回tで微分するとこれになります。逆に...
つまり、バネの例と積分する部分の計算はなにもかわりません。
•DEG→RADの仕方は? ほぼ全員、1回目のレポートで、自分の...
たとえば以下に示すような、各自の工夫を期待します(必須では...
•振り子の長さや初期角度は簡単に変えられるようにする。
•θが小さいときは、sin(√(g/l)t)のただの正弦波で表すことが...
•初期角度θ0と周期の関係を集めて(一種の簡単な実験)、グラフ...
•θ-tのグラフだけではなく、振り子の運動そのものに着目して...
•エネルギー保存則が成り立っているか、計算してみる。
•空気抵抗が働いてもいいかも。
このページここで終了。
&br;
//http://www.mech.tohoku-gakuin.ac.jp/rde/contents/course...
ページ名: