自分のブログログアウトねむねこ幻想郷ネムネコ・ファミリーのために、大昔に書いたファンタジー小説を公開。数学のサイトではありません。[ネコ騙し数学]表計算ソフトで1次元調和振動子を解くためのスプレッドシート1次精度のオイラー法を用いた解法のスプレッドシート(閲覧用)https://docs.google.com/spreadsheets/d/1X8VHH3I08SrWk7EK86SCBOUh3erGy6xJ2fLUvq4G_CQ/pubhtmlスプレッドシートのダウンロード(ここ↓をクリックするとPCにダウンロードされるので注意!!)https://docs.google.com/spreadsheets/d/1X8VHH3I08SrWk7EK86SCBOUh3erGy6xJ2fLUvq4G_CQ/pub?output=odsシンプレクティック法を用いた解法のスプレッドシート(閲覧用)https://docs.google.com/spreadsheets/d/1JDmMdnqOV5peAmmOVgFGZjtWwTu6aoSmiMocmmM9oFM/pubhtmlスプレッドシートのダウンロード(ここ↓をクリックするとPCにダウンロードされるので注意!!)https://docs.google.com/spreadsheets/d/1JDmMdnqOV5peAmmOVgFGZjtWwTu6aoSmiMocmmM9oFM/pub?output=ods閲覧用をクリックすると、例えば、次のようなものが表示されるにゃ。マクロを一切使っていないのでコンピュータウイルスに感染することはないけれど、ネムネコがオレのPCに勝手にスプレッドシートをダウンロードしたと文句をつけられるのは嫌なので、スプレッドシートのダウンロードには注意するにゃ。クリックするとダウンロードのフォルダーにods形式のスプレッドシートがダウンロードされるケロ。これをダウンロードすると、どのように計算しているかが分かると思うにゃ。[ネコ騙し数学]オイラー法とシンプレクティック法を用いた1次元調和振動子の解法1次元調和振動子のハミルトニアン(力学的なエネルギー)はである。ここで、であり、qはバネの自然長からの変位、pは運動量、(1)式のp2/2は運動エネルギー、q2/2はバネの位置エネルギー(ポテンシャルエネルギー)に相当する。(2)式は、と書き換えることが可能で、(2)式とバネの単振動に関するニュートンの運動方程式と同一のものである。時刻t=0における初期条件を(q,p)=(1,0)とすれば、(2)、(3)式の解はである。1次精度のオイラー法を用いて解くには、つまり、と近似し、これを逐次繰り返して計算すればよい。t=0のとき、(q,p)=(1,0)であれば、Δt=0.1とすれば、この計算結果をもとに、といった具合に計算できる。この計算は非常にシンプルな繰り返し計算だから、プログラムを作らずとも、表計算ソフトで簡単に計算が可能で、以下にt=0における初期条件を(q,p)=(1,0)としたものの計算結果を示す。このように非常に粗い近似に基づく計算でありながら、比較的によく計算できていることがわかるだろう。時間の経過とともにqとpの厳密解と1次精度のオイラー法を用いた数値解との食い違いが大きくなるのは、近似計算を繰り返すたびに誤差が伝播し、蓄積するため。このため、オイラー法の誤差は雪だるま式に増加し、計算はやがて発散する。次に、オイラー法を用い、Δt=0.1とし、t=10まで数値計算結果を(q,p)のグラフとして表したものを以下に示す。解いた問題は、力学的エネルギー、つまり、ハミルトニアンHが一定であり、(q,p)は単位円上に存在しなければならないのだが、時間の経過とともに軌道q2+p2=1からの逸脱が大きくなってゆく。オイラー法を用いた近似計算法はエネルギー保存則を満たす解法ではなないということだ。オイラー法のこの欠点を改良したものにシンプレクティック(Sympretic)法と呼ばれるものがある。オイラー法で用いる計算式(4)を次のように改良したもの。この僅かな改良によってエネルギー保存則からの逸脱がかなり改善され、オイラー法の宿命的な欠点である時間の経過に伴う誤差の伝播、蓄積が解消される。以下に、表計算ソフトを用い、シンプレティック法を用いて解いた計算結果を示す。比較参照のために、シンプレクティク法とオイラー法を用いて解いたqと厳密解をのグラフを以下に示す。シンプレティック法による数値解を(q,p)のグラフで表したものを示す。オイラー法とは異なり、シンプレクティック法を用いた数値解はかなり良くq2+p2=1の円周上に存在していることがわかるだろう。シンプレクティク法とオイラー法を用いて解いたq,pを元に計算したハミルトニアンHと厳密な値1/2との相対誤差の時刻による推移をグラフに示す。このグラフが示すように、1次精度のシンプレクティク法は完全にエネルギー保存則を満たすわけでないが、振動をしつつも常にある誤差の範囲内に収まり、オイラー法のように指数関数的に誤差が増大することはない。参考までに2次精度をのシンプレクティク方による計算結果も示してある。数値計算ではよくあることなのだけれど、(4)式を少し変え(5)式にし計算するだけで計算結果が劇的に変化するのだから、数値計算は奥が深いよね。思いませんか?そして、少しでもこうした数値計算に興味を持ってもらえれば幸せです。これらの議論をより正確にするためには、ニュートン力学の進化形である解析力学などの知識が必要になるので、感覚的に理解できれば十分ではないでしょうか。参考までに2次精度のシンプレクティック法で解いた計算結果は以下に示す。合成関数の微分法1定理11関数f(u)が区間Iで微分可能、関数g(x,y)が領域Dで偏微分可能かつならば、合成関数F(x,y)=f(g(x,y))はDで偏微分可能でである。変数yを固定すると、g(x,y)はxのだけの一変数関数φ(x)と考えることができる。1変数関数の合成関数の微分法よりところで、よって、yに関する偏微分の場合も同様。(証明終了)上の定理は、u=g(x,y)、z=f(u)とし、この合成関数z=f(g(x,y))の偏微分をとしたほうが、計算するときにも間違えることはないだろう。次の関数の偏導関数を求めよ。【解】u=x2+y2とおくと、(1)、(2)、(3)の関数はと書き換えられる。次の関係式が成り立つことを示せ。u=ax+byとおくと、z=f(ax+by)=f(u)。u=xyとおくと、z=f(xy)=f(u)。u=x2–y2とおくと、z=f(x2–y2)=f(u)。u=y/xとおくと、z=f(y/x)=f(u)。z=f(x–ct)+g(x+ct)(cは定数)のとき、であることを示せ。v=x+ctとおくと、z=f(u)+g(v)。(解答終了)は波動方程式と呼ばれるもので、z=f(x–ct)+g(x+ct)はその解である。数値積分法を教えてもらったので、バタバタとプログラムを作って計算してみた。数値計算につかったモデルは、調和振動子で、微分方程式はで、t=0のときの初期条件を(q,p)=(1,0)として、1次精度のオイラー法と1次精度のシンプレクティック法を使って、この微分方程式を解いてみた。進めるタイムステップΔt=0.1として数値的に解いた。(1)式のHはハミルトニアンと呼ばれるもので、この場合、運動エネルギーp2/2とバネの位置エネルギーq2/2の和になる。力学的エネルギーが一定だから、(2)を数値的に解いた(q,p)をプロットすると、本来、数値解はの円周上に存在しなければならない。1次精度のオイラー法を用いて解いた数値解は時間の経過とともにp2+q2=1の円から離れていくんだケロよ。エネルギー保存則が満たされていない!!軌道からどんどん離れてゆくから、時間の経過とともにエネルギーが増えていく!!1次精度のオイラー法を用いた数値計算結果これに対し、シンプレクティック法は、多少振動するけれど、力学的なエネルギーが保存されるんだケロ。1次精度のオイラー法だと1次精度のシンプレクティック法だとだけの差なのだけれど、この僅かな差でこれほど計算結果が異なるというのは正直驚きだにゃ。ネムネコは、日々進化しているにゃ。2定点(−a,0)、(a,0)からの距離の積がa2である点Pの軌跡を求めよ。【解】2定点をA(−a,0)、B(a,0)とし、点Pの座標を(x,y)とすると、また、よって、(解答終)問題で求めた点Pの軌跡(曲線)をレムニスケート(lemniscate)という。x=rcosθ、y=rsinθとおき、(1)を極座標を用いて書き直すと、レムニスケートはになる。(2)を用いると、レムニスケートによって囲まれる領域の面積Sをと求めることができる。問題2関係式(x2+y2)2=2a2(x2–y2)(a>0)で定まるxの陰関数yの極値を求めよ。極値を取る点ではy'=0だからx=0では極値を取らないから、x2+y2=a2でなければならない。そこで、(x2+y2)2=2a2(x2–y2)にy2=a2–x2を代入すると、(解答終(?))上の計算からわかるように、極値を取る点はレムニスケート(x2+y2)2=2a2(x2–y2)と原点を中心とする半径aの円x2+y2=a2の交点である。原点(0,0)における接線はy=x、y=−xの2本で、原点は結節点である。おける2次微分係数y''の値を調べる必要があるけれど、計算が複雑になるので、この吟味は省略した。倍角公式よりこれを(3)式に代入すると0≦2θ<π/2=90°とすると、θ=π/6=30°のときに、yは最大値a/2をとる。曲線x2–y2=a2(a>0)の(x0,y0)における接線を求めよ。(1)で求めた接線に原点Oから下ろした垂線の足の軌跡を求めよ。焦点とするレムニスケートである。(解答終)この他にも、レムニスケートは楕円関数などとも深い関係があるけれど、これは「ねこ騙し数学」の現在のレベルを越えてしまうので、この点については触れないことにする。デカルトの正葉線で定まるxの関数yの極値を求めよ。両辺を3で割ると、yが極値を取るときy'=0だから、これを(1)式に代入すると、x=0のとき、y=0(この点(0,0)では(3)式の分母が0になるので微分不可能)。極値を取る点ではy'=0だから、したがって、a>0のときはy''<0で極大、a<0のときはy''>0で極小。取ることだけは確かである。偏微分を使うともう少しスッキリ解くことができるけれど、高校の微分の範囲でデカルトの正葉線で定まる関数の極値を求めることができた。この曲線は、x=0のときy=0になるので、x≠0とし、t=y/xとおくと、(1)式はしたがって、t=0のとき、(x,y)=(0,0)となるので、この曲線はとパラメータ(媒介変数)tを使って表すことができる。パラメータtで表現された(4)を使うと、dy/dx、d2y/dx2が次のように求められるので、これを使って極値や曲線の凹凸を調べることができる。ときに正になるので曲線は下に凸である。参考までに、a>0のときの、xとt、yとtの関係をグラフで示す。曲線をパラメータ表示すると、xとyの直接の対応が失われるので、直観的に理解しづらい!!このグラフを見ると明らかなように、この曲線上の点を(x,y)とすると、だから、|x|→∞のときt→−1である。この曲線の漸近線をy=mx+nとすると、したがって、デカルトの正葉線の漸近線の方程式は、である。原点O(0,0)における接線はy=0、x=0で、原点Oは結節点である。極座標x=rcosθ、y=rsinθを用いると、デカルトの正葉線はr≠0のとき、と表せる。この結果を用いると、ループ、曲線で囲まれた領域の面積Sは次のように求めることができる。Dで偏微分可能のとき、f(x,y)は2回偏微分可能であるといい、fの第2次偏導関数をと定義する。(注)と順序が違うことに注意が必要。次の第2次偏導関数を求めよ。