以下のインプットを実行するには、それをマウスで選び、Shiftを押さえながらReturnを押すか、Enter(キーボードの一番右下）を押して下さい。

単振動

三角関数

まず三角関数の基本的な性質について復習しておこう．

指数関数と三角関数

指数関数と三角関数の間には以下の関係がある(オイラーの公式）．

Exp[I x]==Cos[x]+I Sin[x]

((または，e)^(i x))==cos(x)+i sin(x)

これは、例えば、指数関数、三角関数のTaylor展開をすれば推測できる．即ち、

Series[Exp[I x] , {x, 0, 10}]

Series[Cos[ x] , {x, 0, 10}]

Series[Sin[ x] , {x, 0, 10}]

これらより(e^(i x))==Cos[x]+i Sin[x]であることがわかる．Mathematicaで証明するには、

Simplify[Exp[I x] == Cos[x]+I Sin[x]]

問１  (e^((i π)/2))==i であることを証明しなさい．

加法定理

TrigExpand[ Sin[α+β]]

すなわち

Sin[α+β]==Cos[β] Sin[α]+Cos[α] Sin[β]

ということを意味する．これが真であるかどうかを調べるには

Simplify[Sin[α+β]==Cos[β] Sin[α]+Cos[α] Sin[β]]

となり真実であることがわかる．cosine の場合は、

TrigExpand[ Cos[α+β]]

すなわち

Cos[α+β]==Cos[α] Cos[β]-Sin[α] Sin[β]

である．

これらを一度に証明するためには上の

(e^(i x))==cos(x)+i sin(x)

を持ちいる．即ち、

Exp[I(α+β)] ==Cos[α+β]+I Sin[α+β]

であるが、

TrigExpand[ExpToTrig[Exp[I α] Exp[I β]]]

より、実部と虚部を比べれば、加法定理が正しいことわかる．

問2　以下の式をTrigExpandを使って展開するとどのような結果がえられるか、実際計算して見なさい．
a)  Tan[α+β]
b) Sec[α+β]
c) Csc[α+β]
d) Sin[α+β +γ]
e) Cos[α-2β]

三角関数の複雑な表現を簡単にするためにはSimplify、またはFullSimplifyを使えばよい．例えば

Simplify[Cos[α+2β]+Cos[α-2β]]

Simplify[((Cos[α])^3) ((Cos[β])^2)+3 ((Cos[α])^2) ((Cos[β])^2) Sin[α]-3 Cos[α] ((Cos[β])^2) ((Sin[α])^2)-((Cos[β])^2) ((Sin[α])^3)+2 ((Cos[α])^3) Cos[β] Sin[β]+6 ((Cos[α])^2) Cos[β] Sin[α] Sin[β]-6 Cos[α] Cos[β] ((Sin[α])^2) Sin[β]-2 Cos[β] ((Sin[α])^3) Sin[β]-((Cos[α])^3) ((Sin[β])^2)-3 ((Cos[α])^2) Sin[α] ((Sin[β])^2)+3 Cos[α] ((Sin[α])^2) ((Sin[β])^2)+((Sin[α])^3) ((Sin[β])^2)]

問3　以下の式をSimplifyを使って式を簡略化するとどのような結果がえられるか．
a)  ((Cos[α])^3)-3 Cos[α] ((Sin[α])^2)
b) ((Cos[β])^2)+2 Cos[β] Sin[β]-((Sin[β])^2)
c) 4 ((Cos[β])^3) Sin[β]-4 Cos[β] ((Sin[β])^3)
d) ((4 Cos[β] (Cos[β]-Sin[β]) Sin[β] (Cos[β]+Sin[β]))/((((Cos[β])^2)-2 Cos[β] Sin[β]-((Sin[β])^2)) (((Cos[β])^2)+2 Cos[β] Sin[β]-((Sin[β])^2))))
e) ((Cos[α])^2)+((Sin[α])^2)

単振動

微分方程式の解

単振動を表す運動方程式は

(((d^2)x)/(d(t^2)))=-(ω^2) x

と書ける．（上の表現をマウスで選んでセルメニューの形式変換の中の選択を選べばどのように表現がかわるかやってみよう．）ただし、全微分(Dt)を使うとMathematicaは定数まで変数として解釈してしまい問題が起こるので, 偏微分(D)を以後使う．

単振動の解は正弦関数である．これを

x = A Sin[ω t +δ]

であるとすると

D[x,{t,2}] == -ω^2 x

より

x = A Sin[ω t +δ]

が微分方程式を満足する解であることがわかる．A, δは定数で初期条件によって決まる．

実際に微分方程式を解いてみよう．微分方程式を解くときはDSolveを用いる．

DSolve[D[x[t],{t,2}] == -ω^2 x[t], x[t], t]

アウトプットがおかしいのはxが変数ではなく　上で定義した値を持っているからである．x を入力すると

x

まず、この値を消しておく．

Clear[x]

でも良いし、また、

x=.

でも良い．

x を入力すると

x

となり、x に対して何も値が定義されていないことがわかる．もう一度、計算すると

DSolve[D[x[t],{t,2}] == -ω^2 x[t],x[t], t]

というアウトプットを得る．これは一般解が二つの定数C[1], C[2]を用いてコサインをサインの和として書き表せることを意味する．

TrigExpand[A Sin[δ+t ω]]

であるから

C[1] == A Cos[δ]

C[2]==A Sin[δ]

であることがわかる．もちろん、二つの定数C[1], C[2]は初期条件によって決まる．

また、微分方程式

(((d^2)x)/(d(t^2)))==-(ω^2) x

は

(((d^2)x)/(d(t^2)))==((dv)/(dt))==((dv)/(dx)) ((dx)/(dt))==v ((dv)/(dx))

であるから

v ((dv)/(dx))==-(ω^2) x

とも書ける．これを積分すると

((v^2)/2)+(((ω^2) (x^2))/2)==(Energy/m)

となる．ここでEnergy, m は定数である（それぞれエネルギー（運動エネルギーと位置エネルギーの和）、質量に対応している）．以下、少し分かりにくいかもしれないが、その計算の手順を示す．

Solve[((v^2)/2)+(((ω^2) (x^2))/2)==(Energy/m), v]

従って、

(1/(Sqrt[((2 Energy)/m)-(x^2) (ω^2)]))dx == dt

であるから左辺を積分しよう．

∫(1/(Sqrt[((2 Energy)/m)-(x^2) (ω^2)]))\[DifferentialD]x

この結果はすこし複雑に見える．FullSimplifyして

FullSimplify[%]

となり、わかりやすい結果を得る．次 x について解いてみよう．ここで積分定数をconstant/ωとする．

Solve[((ArcTan[((x ω)/(Sqrt[((2 Energy)/m)-(x^2) (ω^2)]))])/ω) == t+ constant/ω, x]

まだ複雑に見えるので、FullSimplifyを使う．

FullSimplify[%]

もう少し簡単にする．

FullSimplify[%, ω >0&& Sec[constant+t ω]>0]

となり、求めたい結果を得る．変位 x 、速度v は

x=x/.First[%]

v = D[x, t]

全エネルギーを計算すると

((m (v^2))/2)+(1/2) m (ω^2) (x^2)

これを簡単にすると

Simplify[%]

となり解が正しいことがわかる．

図を描く．

ここで位置、速度、運動エネルギー、位置エネルギー等を時間の関数としてプロットしてみよう．いま、位置　ｘは

x

であるが、簡単のために

ω=1;
Energy=1;
m=1;
constant=0;

と仮定する．
まず、位置 x を時間の関数としてプロットしてみよう．

Plot[x, {t, 0,10}]

Optionを付けると、

Plot[x, {t, 0,10},AxesLabel->{時間 , 変位}, PlotStyle->Hue[1]]

同様に、速度については

Plot[v, {t, 0,10},AxesLabel->{時間 ,  速度}, PlotStyle->Hue[.8]]

速度は変位に比して位相が９０度遅れているのが分かる．加速度は

Plot[Evaluate[D[v,t]], {t, 0,10},AxesLabel->{時間 ,  加速度}, PlotStyle->Hue[.4]]

はさらに位相が９０度遅れているのが分かる．

さて、運動エネルギーはどうであろうか．

graphKinetic=Plot[(1/2) m (v^2),{t,0,10},AxesLabel->{時間,運動エネルギー},PlotStyle->Hue[0.9]]

ここでアウトプットのグラフをgraphKineticと呼んだことに注意．位置エネルギーは

graphPotential=Plot[(1/2) m (ω^2) (x^2),{t,0,10},AxesLabel->{時間,位置エネルギー},PlotStyle->Hue[0.7]]

全エネルギーは

graphTotalEnergy=Plot[Evaluate[(1/2) m (ω^2) (x^2)+(1/2) m (v^2)],{t,0,10},AxesLabel->{時間,エネルギー},PlotStyle->Hue[0.5], PlotRange->{1-10^(-15),1+10^(-15)}]

以上の図をまとめて描くと、

Show[graphTotalEnergy,graphKinetic,graphPotential, PlotRange->{-0.1, 1.1}]

問４　上の図の物理的な解釈を述べよ．(運動エネルギーと位置エネルギー，力学的エネルギーの関係は？）

ここで教科書で説明してある相空間図を描いておく．

ParametricPlot[{x, v}, {t, 0, 10},AxesLabel->{変位,速度},PlotStyle->Hue[0.95], Ticks->None]

いろいろなエネルギーの相空間図を描きたいときは、

Energy=.

ParametricPlot[Evaluate[Table[{x, v},{ Energy,.01,1,.05}]], {t, 0, 10},AxesLabel->{変位,速度},PlotStyle->Table[Hue[Energy],{ Energy,.1,1,.1}], Ticks->None]

問５　教科書の課題１−１，１−２を解け．

問６　(optional) 教科書の課題１−３，１−４を解け．