表計算ソフトで１次元調和振動子を解くためのスプレッドシート [ネコ騙し数学]
https://docs.google.com/spreadsheets/d/1X8VHH3I08SrWk7EK86SCBOUh3erGy6xJ2fLUvq4G_CQ/pubhtml
https://docs.google.com/spreadsheets/d/1X8VHH3I08SrWk7EK86SCBOUh3erGy6xJ2fLUvq4G_CQ/pub?output=ods
https://docs.google.com/spreadsheets/d/1JDmMdnqOV5peAmmOVgFGZjtWwTu6aoSmiMocmmM9oFM/pubhtml
https://docs.google.com/spreadsheets/d/1JDmMdnqOV5peAmmOVgFGZjtWwTu6aoSmiMocmmM9oFM/pub?output=ods
ネムネコがオレのPCに勝手にスプレッドシートをダウンロードした
と文句をつけられるのは嫌なので、スプレッドシートのダウンロードには注意するにゃ。
クリックするとダウンロードのフォルダーにods形式のスプレッドシートがダウンロードされるケロ。これをダウンロードすると、どのように計算しているかが分かると思うにゃ。
オイラー法とシンプレクティック法を用いた１次元調和振動子の解法 [ネコ騙し数学]
オイラー法とシンプレクティック法を用いた１次元調和振動子の解法
１次元調和振動子のハミルトニアン（力学的なエネルギー）は
である。
ここで、
であり、qはバネの自然長からの変位、pは運動量、（１）式のp²/2は運動エネルギー、q²/2はバネの位置エネルギー（ポテンシャルエネルギー）に相当する。
（２）式は、
と書き換えることが可能で、（２）式とバネの単振動に関するニュートンの運動方程式と同一のものである。
時刻t=0における初期条件を(q,p)=(1,0)とすれば、（２）、（３）式の解は
である。
１次精度のオイラー法を用いて解くには、
つまり、
と近似し、これを逐次繰り返して計算すればよい。
たとえば、t=0のとき、(q,p)=(1,0)であれば、Δt=0.1とすれば、
この計算結果をもとに、
といった具合に計算できる。
この計算は非常にシンプルな繰り返し計算だから、プログラムを作らずとも、表計算ソフトで簡単に計算が可能で、以下にt=0における初期条件を(q,p)=(1,0)としたものの計算結果を示す。
このように非常に粗い近似に基づく計算でありながら、比較的によく計算できていることがわかるだろう。時間の経過とともにqとpの厳密解と１次精度のオイラー法を用いた数値解との食い違いが大きくなるのは、近似計算を繰り返すたびに誤差が伝播し、蓄積するため。このため、オイラー法の誤差は雪だるま式に増加し、計算はやがて発散する。
次に、オイラー法を用い、Δt=0.1とし、t=10まで数値計算結果を(q,p)のグラフとして表したものを以下に示す。
解いた問題は、力学的エネルギー、つまり、ハミルトニアンHが一定であり、
(q,p)は単位円上に存在しなければならないのだが、時間の経過とともに軌道q²+p²=1からの逸脱が大きくなってゆく。
つまり、オイラー法を用いた近似計算法はエネルギー保存則を満たす解法ではなないということだ。
オイラー法のこの欠点を改良したものにシンプレクティック（Sympretic）法と呼ばれるものがある。
オイラー法で用いる計算式（４）を次のように改良したもの。
この僅かな改良によってエネルギー保存則からの逸脱がかなり改善され、オイラー法の宿命的な欠点である時間の経過に伴う誤差の伝播、蓄積が解消される。
以下に、表計算ソフトを用い、シンプレティック法を用いて解いた計算結果を示す。
比較参照のために、シンプレクティク法とオイラー法を用いて解いたqと厳密解をのグラフを以下に示す。
さらに、シンプレティック法による数値解を(q,p)のグラフで表したものを示す。
オイラー法とは異なり、シンプレクティック法を用いた数値解はかなり良くq²+p²=1の円周上に存在していることがわかるだろう。
シンプレクティク法とオイラー法を用いて解いたq,pを元に計算したハミルトニアンHと厳密な値1/2との相対誤差
の時刻による推移をグラフに示す。
このグラフが示すように、１次精度のシンプレクティク法は完全にエネルギー保存則を満たすわけでないが、振動をしつつも常にある誤差の範囲内に収まり、オイラー法のように指数関数的に誤差が増大することはない。
参考までに２次精度をのシンプレクティク方による計算結果も示してある。
数値計算ではよくあることなのだけれど、（４）式を少し変え（５）式にし計算するだけで計算結果が劇的に変化するのだから、数値計算は奥が深いよね。
そう思いませんか？
そして、少しでもこうした数値計算に興味を持ってもらえれば幸せです。
なお、これらの議論をより正確にするためには、ニュートン力学の進化形である解析力学などの知識が必要になるので、感覚的に理解できれば十分ではないでしょうか。
参考までに２次精度のシンプレクティック法で解いた計算結果は以下に示す。
第１０回 合成関数の微分法１ [ネコ騙し数学]
第１０回 合成関数の微分法１
定理１１
関数f(u)が区間Iで微分可能、関数g(x,y)が領域Dで偏微分可能かつ
ならば、合成関数F(x,y)=f(g(x,y))はDで偏微分可能で
である。
【証明】
変数yを固定すると、g(x,y)はxのだけの一変数関数φ(x)と考えることができる。
１変数関数の合成関数の微分法より
ところで、
よって、
yに関する偏微分の場合も同様。
（証明終了）
上の定理は、u=g(x,y)、z=f(u)とし、この合成関数z=f(g(x,y))の偏微分を
としたほうが、計算するときにも間違えることはないだろう。
問１ 次の関数の偏導関数を求めよ。
【解】
u=x²+y²とおくと、（１）、（２）、（３）の関数は
と書き換えられる。
z = f(u)とおき、uで微分すると
u=x²+y²の偏導関数は
だから、
（解答終了）
問２ 次の関係式が成り立つことを示せ。
（１） z=f(ax+by)のとき、
（２） z=f(xy)のとき
（３） z=f(x²–y²)のとき
（４） 
のとき
【解】
（１） u=ax+byとおくと、z=f(ax+by)=f(u)。
また、
。
よって、
（２） u=xyとおくと、z=f(xy)=f(u)。
また、
。
（３） u=x²–y²とおくと、z=f(x²–y²)=f(u)。
また、
だから、
（４） u=y/xとおくと、z=f(y/x)=f(u)。
また、
だから、
（解答終）
問３ z=f(x–ct)+g(x+ct)（cは定数）のとき、
であることを示せ。
【解】
u=x–ct 、v=x+ctとおくと、z=f(u)+g(v)。
また、
だから、
よって、
である。
（解答終了）
は波動方程式と呼ばれるもので、z=f(x–ct)+g(x+ct)はその解である。
ネムネコ、シンプレクティック法を使って微分方程式の数値積分をする！！ [ネコ騙し数学]
かつてのネムネコ・ファミリー――某Q&Aサイト上のみでネムネコと深い交流のあったヒトたち――であるddt³さんから、シンプレクティック（Symplectic）法と呼ばれる数値計算法、微分方程式の数値積分法を教えてもらったので、バタバタとプログラムを作って計算してみた。
数値計算につかったモデルは、調和振動子
で、微分方程式は
で、t=0のときの初期条件を(q,p)=(1,0)として、１次精度のオイラー法と１次精度のシンプレクティック法を使って、この微分方程式を解いてみた。
進めるタイムステップΔt=0.1として数値的に解いた。
（１）式のHはハミルトニアンと呼ばれるもので、この場合、運動エネルギーp²/2とバネの位置エネルギーq²/2の和になる。力学的エネルギーが一定だから、（２）を数値的に解いた(q,p)をプロットすると、本来、数値解は
の円周上に存在しなければならない。
しかし、１次精度のオイラー法を用いて解いた数値解は時間の経過とともにp²+q²=1の円から離れていくんだケロよ。つまり、エネルギー保存則が満たされていない！！
軌道からどんどん離れてゆくから、時間の経過とともにエネルギーが増えていく！！
１次精度のオイラー法を用いた数値計算結果
これに対し、シンプレクティック法は、多少振動するけれど、力学的なエネルギーが保存されるんだケロ。
１次精度のシンプレクティック法を用いた計算結果
数値積分法は、
１次精度のオイラー法だと
１次精度のシンプレクティック法だと
だけの差なのだけれど、この僅かな差でこれほど計算結果が異なるというのは正直驚きだにゃ。
レムニスケート [ネコ騙し数学]
レムニスケート
問題１ ２定点（−a,0）、(a,0)からの距離の積がa²である点Pの軌跡を求めよ。ただし、(a>0)とする。
【解】
２定点をA（−a,0）、B(a,0)とし、点Pの座標を(x,y)とすると、
また、
よって、
（解答終）
問題で求めた点Pの軌跡（曲線）
をレムニスケート（lemniscate）という。
x=rcosθ、y=rsinθとおき、（１）を極座標を用いて書き直すと、レムニスケートは
になる。
さらに、（２）を用いると、レムニスケートによって囲まれる領域の面積Sを
と求めることができる。
問題２
関係式(x²+y²)²=2a²(x²–y²)（a>0）で定まるxの陰関数yの極値を求めよ。
【解】
(x²+y²)²=2a²(x²–y²)の両辺をxで微分すると、
極値を取る点ではy'=0だから
x=0では極値を取らないから、x²+y²=a²でなければならない。
そこで、(x²+y²)²=2a²(x²–y²)にy²=a²–x²を代入すると、
（解答終（？））
上の計算からわかるように、極値を取る点はレムニスケート(x²+y²)²=2a²(x²–y²)と原点を中心とする半径aの円x²+y²=a²の交点である。
原点(0,0)における接線はy=x、y=−xの２本で、原点は結節点である。
うるさいことを言えば
で極値になることを示すために、における２次微分係数y''の値を調べる必要があるけれど、計算が複雑になるので、この吟味は省略した。
ただ、微分積分を使わずに、極座標を用いて次のように解くこともできるだろう。
【別解】
三角関数の倍角公式より
これを（３）式に代入すると
0≦2θ<π/2=90°とすると、θ=π/6=30°のときに、yは最大値a/2をとる。
このとき、
だから、
この曲線はx軸、y軸、原点に関して対称だから、極値は
（別解終（？））
問題３
（１） 曲線x²–y²=a²（a>0)の(x₀,y₀)における接線を求めよ。
（２） （１）で求めた接線に原点Oから下ろした垂線の足の軌跡を求めよ。
【解】
（１） x²–y²=a²の両辺をxで微分すると、
よって、曲線上の点を(x₀,y₀)とするとき、y₀≠0のとき、接線の方程式は
接線の方程式は、y=0、x=−aのときx=−a、y=0、x=aのときx=aとなるが、この場合も⑨式で表せるので、
（２） 曲線上の点(x₀,y₀)における曲線の接線の方程式は、x₀x–y₀y=a²だから、原点Oを通りこの接線に垂直な直線の方程式は、
したがって、垂線の足の座標は連立方程式
の解で
よって、
曲線x²–y²=a²上の点(x₀,y₀)をx₀=ρcosθ、y₀=ρsinθとすると、
これを⑨³に代入すると、
となり、軌跡は
を焦点とするレムニスケートである。
（解答終）
この他にも、レムニスケートは楕円関数などとも深い関係があるけれど、これは「ねこ騙し数学」の現在のレベルを越えてしまうので、この点については触れないことにする。
デカルトの正葉線 [ネコ騙し数学]
デカルトの正葉線
問題 デカルトの正葉線
で定まるxの関数yの極値を求めよ。
【解】
（１）式の両辺をxで微分すると、
両辺を３で割ると、
yが極値を取るときy'=0だから、
これを（１）式に代入すると、
x=0のとき、y=0（この点(0,0)では(3)式の分母が0になるので微分不可能）。
のとき、
極値の判定をするために、（２）式をxで微分しy''を求めると、
極値を取る点ではy'=0だから、
したがって、a>0のときはy''<0で極大、a<0のときはy''>0で極小。
（解答終）
(0,0)を極値と認めるかどうかという微妙な問題があるけれど、x³+y³–3axy=0で定まるxの関数（陰関数)yがのときに極値
を取ることだけは確かである。
また、x³+y³–3axy=0はxとyを入れ替えても成立するので――この曲線はy=xに関して対称――,x³+y³–3axy=0で定まるyの関数xが極値を取る点は曲線上の点
となる。
偏微分を使うともう少しスッキリ解くことができるけれど、高校の微分の範囲でデカルトの正葉線で定まる関数の極値を求めることができた。
この曲線は、x=0のときy=0になるので、x≠0とし、t=y/xとおくと、（１）式は
したがって、
t=0のとき、(x,y)=(0,0)となるので、この曲線は
とパラメータ（媒介変数）tを使って表すことができる。
パラメータtで表現された（４）を使うと、dy/dx、d²y/dx²が次のように求められるので、これを使って極値や曲線の凹凸を調べることができる。
a>0のとき、d²y/dx²は、
で負となるので曲線は上に凸、
のときに正になるので曲線は下に凸である。
参考までに、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は次のように求めることができる。
ここで、tanθ=tとおくと
第９回 高階偏導関数 [ネコ騙し数学]
第９回 高階偏導関数
領域D上で定義された偏微分可能な関数f(x,y)の導関数
がDで偏微分可能のとき、f(x,y)は２回偏微分可能であるといい、fの第２次偏導関数を
と定義する。
（注）
と順序が違うことに注意が必要。
問１ 次の第２次偏導関数を求めよ。
【解】
（１） z=f(x,y)=x³+xy²+y⁴とおくと、
だから、
（２） z=f(x,y)=sin(xy)とおくと、
だから、
（解答終了）
領域D上で定義された関数f(x,y)が第n次導関数を持ち、それらがDで連続であるとき、fはn回連続微分可能または
級であるという。任意のn∈Nに対して級である関数を無限回連続微分可能または
級という。
定理
関数f(x,y)がC²級ならば、
が成り立つ。
【証明】
とおき、平均値の定理を用いると
同様に、
従って、
fはC²級だから２次偏導関数は連続で、(h,k)→(0,0)のとき、
になる。
（証明終了）
問１の（１）、（２）ともにになっているのは偶然ではなく、f(x,y)がC²級だからである。また、C²級ならば、が成立するので、片方だけを計算すればよい。
問題１ ２回偏微分可能な関数f(x,y)に
を対応させる写像を
であらわし、ラプラス微分作用素またはラプラシアンという。Δf=0をラプラスの微分方程式、その解を調和関数という。
次の関数z=f(x,y)が調和関数であるかどうか調べよ。
【解】
（１） 
より
よって、調和関数である。
（２） 
より
だから、調和関数ではない。
（３） 
より
よって、調和関数である。
（解答終了）
問題２ f(x,y)をR²の開区間(a,b)×(c,d)で定義された関数とする。次のことを証明せよ。
（１） fがxで偏微分可能で
ならば、fはyだけの関数である。
（２） 
ならば、fは定数である。
（３） 
かつ
が連続ならば、fはxだけの関数とyだけの関数の和である。
【証明】
（１） yを固定すると、xだけの関数g(x)=f(x,y)は(a,b)上で
だから定数である。この定数をφ(y)とおくと、f(x,y)=φ(y)である。
（２） だから（１）よりf(x,y)=φ(y)。
よって、
したがって、φ(y)は定数である。
（３） 
だから、
はxだけの関数h(x)に等しい。
x₀∈(a,b)の定点とする。
はxについて連続だから
ここで、
とおけば、
（証明終了）
（注）
(a,b)×(c,d)は開区間(a,b)と開区間(c,d)の直積のことで
である。
また、は、Dのすべての点(x,y)で
の意味で、「fが定数」とはDのすべての点(x,y)でf(x,y)=定数、つまり、関数fが定数関数であることを表す。
「関数f(x,y)=定数」ではなく、「関数f=定数」と書く方が正式！！
完全微分形微分方程式２ [ネコ騙し数学]
完全微分形微分方程式２
定義
１階の全微分方程式
の左辺がある関数f(x,y)の全微分、すなわち、
であるとき、全微分方程式を完全形であるという。
定理
関数P(x,y)、Q(x,y)の偏導関数が領域
で連続であるとき、次の条件は同値である。
（１） P(x,y)dx+Q(x,y)dy=0は完全形である。
（２）
【証明】
（１）⇒（２）
P(x,y)dx+Q(x,y)dy=0は完全形だから、
になる関数f(x,y)が存在する。
よって、
仮定より
は連続だから
も連続で
。
したがって、
（２）⇒（１）
点(x₀,y₀)を定義域内の点とする。
だから、
これを
に代入すると、
Pの偏導関数は連続だから、微分と積分の順序の交換が可能。
よって、
そこで、
と定めると、
したがって、
になり、P(x,y)dx+Q(x,y)dy=0は完全形である。
（証明終）
以上のことより、完全型微分方程式
の一般解z=f(x,y)は
である。
問 次の微分方程式の一般解を求めよ。
【解】
（１） P(x,y)=x、Q(x,y)=yとおくと、
だから、この微分方程式は完全形。
x₀=0、y₀=0とすると、公式より
（２） P(x,y)=x²–y 、Q(x,y)=y²–xとおくと、
よって、この微分方程式は完全形。
x₀=0、y₀=0とすると、公式より
（解答終）
公式を用いれば、このように機械的に計算することができる。
本によっては、完全形微分方程式の一般解z=f(x,y)を
としているものがあるので、これを証明することにする。
【略証】
これを、
に代入すると、Pの偏導関数は連続だから微分と積分の順序の交換が可能だから、
したがって、
よって、一般解は
（略証終）
確かに、一般解はこのようになるけれど、積分を３回、さらに偏微分を１回しなければいけないので、この公式から一般解を求めるのは大変。だから、使わないほうが身のためだケロ！！
そもそも、前回紹介した完全形の全微分方程式の解法は、これと同じことをやっている。
したがって、この公式は覚える必要は全く無いように思う。
完全形微分方程式 [ネコ騙し数学]
完全形微分方程式
関数z=f(x,y)の全微分は
である。
問題１ ある関数z=f(x,y)の全微分が（１）、（２）で与えられているとき、関数zを求めなさい。
【解】
（１） zの全微分は
だから、
①をxで積分すると、
これを②に代入すると、
したがって、
（２）
yを固定して、③式をxで積分すると、
これを④式に代入すると、
よって、
（解答終了）
（※） φ(y)はxを含まないyだけの関数。
また、
となるから、
である。
この問題は（全）微分の性質を用いて、次のように解いてもよい。
【別解】
（別解終）
定義 １階の全微分方程式
の左辺がある関数z=f(x,y)の全微分dzに等しいとき、すなわち、
であるとき、完全形であるという。
問題２
関数f(x,y)の全微分が、
であるとき、関係式f(x,y)=C（Cは定数）によって、微分方程式
の解が得られる。
（１） 
ならば
であることを示せ。
（２） 次の方程式が
を満たすことを確かめて、これを解け。
【解】
（１）
が連続であるときだから、
である。
（２a) P=x²–y 、Q=y²–xだから、
だから、yを固定してxで積分すると、
また、
だから、（１）をこれに代入すると、
よって、
(2b) P=2xy、Q=x²–y²とおくと、
だから
にこれを代入すると、
よって、
第８回 全微分に関する問題 [ネコ騙し数学]
第８回 全微分に関する問題
問題１ 次の問に答えよ。
（１） 
の全微分を求めよ。
（２） x、yの関数f、gに関して次のことが成り立つことを示せ。
【解】
（１）
（２）
（解答終了）
問題２ 次の問に答えよ。
（１） ２辺がa、bである長方形の面積をSとする。a、bがそれぞれΔa、Δb変化したときの面積の変化量をΔSとすると、ΔS=bΔa+aΔbであることを示せ。
（２） △ABCが定円に内接しながら微小変化したとき、３辺a、b、cの変化量をΔa、Δb、Δcとすると
が成り立つことを示せ。
【解】
（１） S=abの全微分は
微小な変化だからΔS≒dS。
よって、ΔS=bΔa+aΔbである。
（別解）
ΔaΔbは、bΔaとaΔbより高次の微小項だから、｜Δa｜≪1、｜Δb｜≪1のとき無視することができる。
したがって、ΔS=bΔa+aΔb
（２） 定円の半径をRとすると、正弦定理より
したがって、
同様に、
よって、
da=Δa、db=Δb、dc=Δcだから
（解答終了）
問題２のように、全微分は誤差の評価に応用が可能な場合がある。
問題３ 次の関数は原点(0,0)で微分可能か。
【解】
（１） (x,y)=(0,0)における偏微分係数を求めると、
だから、
(h,k)→(0,0)のとき、
であれば、f(x,y)は(0,0)で微分可能、そうでなければ(0,0)で微分可能でない。
h=t、y=tとき、t→0として(h,k)を(0,0)に近づけると、
したがって、f(x,y)は(0,0)で微分可能ではない。
（２）
(x,y)=(0,0)における偏微分係数を求めると、
だから、
とおくと、｜h｜≦r、｜k｜≦rだから、
したがって、f(x,y)は(0,0)で微分可能である。
（解答終了）
ちなみに、問題３の関数f(x,y)はいずれも(0,0)で微分可能である。
また、（１）のf(x,y)の、(0,0)におけるθ方向の方向微分係数を求めると、
原点における方向微分係数がθの値によって変化し、１つの値に定まらないので、f(x,y)は(0,0)で微分可能ではない。
同様に、（２）の場合の方向微分係数を求めると、
となり、θの値にかかわらず１つの値に定まるので、f(x,y)は(0,0)で微分可能ということになる。
