033_HarmonicOscillators.mw

連性調和振動子

イントロダクション

このワークシートでは、複雑な問題を解き、結果をグラフィカルに見せることによりMapleVのパワーをデモンストレーションします。ワークシートの中のアニメーションはシステムの多様性を表現しています。この例で使われているシステムはバネの振動システムです。

[Inserted Image]

> restart;

システムの数学モデル

> eq1:=m1*diff(x1(t),t$2)+(k1+k12)*x1(t)-k12*x2(t)=0;

eq1 := m1*(diff(x1(t), `$`(t, 2)))+(k1+k12)*x1(t)-k12*x2(t) = 0

> eq2:=m2*diff(x2(t),t$2)-k12*x1(t)+(k12+k2)*x2(t)=0;

eq2 := m2*(diff(x2(t), `$`(t, 2)))-k12*x1(t)+(k12+k2)*x2(t) = 0

> ics:=x1(0)=X1,x2(0)=X2,D(x1)(0)=v1,D(x2)(0)=v2;

ics := x1(0) = X1, x2(0) = X2, D(x1)(0) = v1, D(x2)(0) = v2

システムの解析

> ans:=dsolve({eq1,eq2,ics},{x1(t),x2(t)}):

ブロックのそれぞれについて解を代入します。

> block1:=subs(ans,x1(t)):
block2:=subs(ans,x2(t)):

パラメータを代入します。

m1 = m2
k1 = k12 = k2

x1(0) = 0

x2(0) = 0

v1(0) = 0

v2(0) = v

> assume(k,positive); assume(m,positive);

> b1:=simplify(simplify(subs([k1=k,k2=k,k12=k,m1=m,m2=m,X1=0,X2=0,v1=0,v2=v],block1)));

b1 := 1/6*m^(1/2)*v*(3*sin(k^(1/2)*t/m^(1/2))-3^(1/2)*sin(3^(1/2)*k^(1/2)*t/m^(1/2)))/k^(1/2)

> b2:=simplify(simplify(subs([k1=k,k2=k,k12=k,m1=m,m2=m,X1=0,X2=0,v1=0,v2=v],block2)));

b2 := 1/6*m^(1/2)*v*(3*sin(k^(1/2)*t/m^(1/2))+3^(1/2)*sin(3^(1/2)*k^(1/2)*t/m^(1/2)))/k^(1/2)

残りのパラメータに以下の値を代入します。

m = 1
k = 1

v = 1

> p1:=unapply(subs([m=1,k=1,v=1],b1),t);

p1 := proc (t) options operator, arrow; 1/2*sin(t)-1/6*3^(1/2)*sin(3^(1/2)*t) end proc

> p2:=unapply(subs([m=1,k=1,v=1],b2),t);

p2 := proc (t) options operator, arrow; 1/2*sin(t)+1/6*3^(1/2)*sin(3^(1/2)*t) end proc

> plot([p1(t),p2(t)],t=0..20,color=[red,blue]);

[Plot]

sin(t) および sin(sqrt(3)*t) の周期は決して一致しないので繰り返しは決して起こらないことに注意してください。

バネを描くプロシージャ

> springplot := proc(a,b,n)
 local l,i,p;
 l := (b-a)/n;
 p[0] := plot([[a,0],[a+l,0]],color=black,thickness=2);
 for i from 1 to n-2 do
   p[i] := plot([[a+i*l,0],[a+i*l+l/3,1/2],[a+i*l+2*l/3,-1/2],[a+i*l+l,0]],color=black,thickness=2);
 end do;
 p[n-1] := plot([[a+(n-1)*l,0],[a+(n-1)*l+l,0]], color=black, thickness=2);
 p := plots[display](seq(p[i],i=0..n-1));
end proc:

アニメーション

> with(plots):

Warning, the name changecoords has been redefined

アニメーションを作成するプロシージャ

> oscillator := proc(p1,p2,n::posint,f)
 local end1,end2,a,b,c1,c2,spring1,spring2,spring3,masses,i;
 end1 := plot([[0,-1],[0,1]],color=black,thickness=2);
 end2 := plot([[17,-1],[17,1]],color=black,thickness=2);
 for i from 0 to n-1 do
   a := p1(i*f)+5; b:=p2(i*f)+11;
   c1[i] := polygonplot([[a-1,-1],[a-1,1],[a+1,1],[a+1,-1]],color=red);
   c2[i] := polygonplot([[b-1,-1],[b-1,1],[b+1,1],[b+1,-1]],color=blue);
   spring1[i] := springplot(0,a-1,10);
   spring2[i] := springplot(a+1,b-1,10);
   spring3[i] := springplot(b+1,17,10);
   masses[i] := display({c1[i],c2[i],spring1[i],spring2[i],spring3[i],end1,end2});
 end do:
 display([seq(masses[i],i=0..n-1)],insequence=true,scaling=constrained,axes=none, title=`Spring Oscillating System`);
end proc:

> oscillator(p1,p2,40,1/2);

[Plot]

方程式を手で解きます。

> eq3:=subs([m1=1,m2=1,k1=k,k12=k,k2=k],eq1);

eq3 := (diff(x1(t), `$`(t, 2)))+2*k*x1(t)-k*x2(t) = 0

> eq4:=subs([m1=1,m2=1,k1=k,k12=k,k2=k],eq2);

eq4 := (diff(x2(t), `$`(t, 2)))-k*x1(t)+2*k*x2(t) = 0

システムは線形なので次の形式を仮定します。x1(t) = a*exp(I*(omega*t+phi)) , x2(t) = b*exp(I*(omega*t+phi))

> eq5:=x1(t)=a*exp(I*(omega*t+phi));
eq6:=x2(t)=b*exp(I*(omega*t+phi));

eq5 := x1(t) = a*exp(I*(omega*t+phi))

eq6 := x2(t) = b*exp(I*(omega*t+phi))

仮定の解を連立方程式に代入します。

> eq7:=simplify(subs([eq5,eq6],eq3));
eq8:=simplify(subs([eq5,eq6],eq4));

eq7 := -exp(I*(omega*t+phi))*(a*omega^2-2*k*a+k*b) = 0

eq8 := exp(I*(omega*t+phi))*(-b*omega^2-k*a+2*k*b) = 0

exp(I*(omega*t+phi)) で方程式をわります。

> eq9:=simplify(eq7/exp(I*(omega*t+phi)));
eq10:=simplify(eq8/exp(I*(omega*t+phi)));

eq9 := -a*omega^2+2*k*a-k*b = 0

eq10 := -b*omega^2-k*a+2*k*b = 0

見やすくするために、方程式を再度整理します。

> eq9:=collect(eq9,[a,b]);
eq10:=collect(eq10,[a,b]);

eq9 := (-omega^2+2*k)*a-k*b = 0

eq10 := -k*a+(-omega^2+2*k)*b = 0

この方程式(eq9, eq10)は斉次なので、自明でない解が存在するのには次の行列の行列式がゼロであるという条件が必要です。

> A:=matrix([[coeff(lhs(eq9),a), coeff(lhs(eq9),b)], [coeff(lhs(eq10),a), coeff(lhs(eq10),b)]]);

A := matrix([[-omega^2+2*k, -k], [-k, -omega^2+2*k]])

> with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

Aの行列式を計算します

> det(A);

omega^4-4*omega^2*k+3*k^2

この時点でdet(A)=0を解くこともできますが、この場合は det(A) を上手く因数分解できるので、 それを利用します。

> detA:=factor(det(A));

detA := (3*k-omega^2)*(k-omega^2)

omega を直接解く代わりに omega^2 について解きます。

> omegasquared:=solve(detA=0,omega^2);

omegasquared := k, 3*k

> omega1:=omega^2=omegasquared[1]; omega2:=omega^2=omegasquared[2];

omega1 := omega^2 = k

omega2 := omega^2 = 3*k

omega^2 の値を eq9 および eq10 に代入します。その式を a について解き、b を使った形式で表します。

> subs(omega1,eq9);
subs(omega1,eq10);

k*a-k*b = 0

-k*a+k*b = 0

> solve(%,a);

b

> subs(omega2,eq9);
subs(omega2,eq10);

-k*a-k*b = 0

-k*a-k*b = 0

> solve(%,a);

-b

a, bおよび omega について解った値を eq5 および eq6の中に代入します。

> eq11:=convert(subs([a=a1,b=-a1,omega=sqrt(3*k)],eq5),trig);
eq12:=convert(subs([a=a1,b=-a1,omega=sqrt(3*k)],eq6),trig);

eq11 := x1(t) = a1*(cos(3^(1/2)*k^(1/2)*t+phi)+I*sin(3^(1/2)*k^(1/2)*t+phi))

eq12 := x2(t) = -a1*(cos(3^(1/2)*k^(1/2)*t+phi)+I*sin(3^(1/2)*k^(1/2)*t+phi))

> eq13:=convert(subs([a=a2,b=a2,omega=sqrt(k)],eq5),trig);
eq14:=convert(subs([a=a2,b=a2,omega=sqrt(k)],eq6),trig);

eq13 := x1(t) = a2*(cos(k^(1/2)*t+phi)+I*sin(k^(1/2)*t+phi))

eq14 := x2(t) = a2*(cos(k^(1/2)*t+phi)+I*sin(k^(1/2)*t+phi))

解の実部を取り出します。

> eq15:=lhs(eq11)=evalc(Re(rhs(eq11))+Re(rhs(eq13)));

eq15 := x1(t) = a1*cos(3^(1/2)*k^(1/2)*t+phi)+a2*cos(k^(1/2)*t+phi)

> eq16:=lhs(eq12)=evalc(Re(rhs(eq12))+Re(rhs(eq14)));

eq16 := x2(t) = -a1*cos(3^(1/2)*k^(1/2)*t+phi)+a2*cos(k^(1/2)*t+phi)

速度を計算します (これは初期条件を指定するために必要です)。

> eq17:=v1(t)=diff(rhs(eq15),t);

eq17 := v1(t) = -a1*sin(3^(1/2)*k^(1/2)*t+phi)*3^(1/2)*k^(1/2)-a2*sin(k^(1/2)*t+phi)*k^(1/2)

> eq18:=v2(t)=diff(rhs(eq16),t);

eq18 := v2(t) = a1*sin(3^(1/2)*k^(1/2)*t+phi)*3^(1/2)*k^(1/2)-a2*sin(k^(1/2)*t+phi)*k^(1/2)

この初期条件を使って未知のパラメータを解きます。

> parms :=
        solve(
                {subs(t=0,rhs(eq15))=0,subs(t=0,rhs(eq17))=0, subs(t=0,rhs(eq16))=0,subs(t=0,rhs(eq18))=v},{a1,a2,phi}
        );

parms := {a1 = 1/6*v*3^(1/2)/k^(1/2), phi = 1/2*Pi, a2 = -1/2*v/k^(1/2)}, {phi = -1/2*Pi, a1 = -1/6*v*3^(1/2)/k^(1/2), a2 = 1/2*v/k^(1/2)}

最終的な結果を得るためにパラメータの値を代入します。

> simplify(subs(parms[1],eq15)); simplify(subs(parms[1],eq16));

x1(t) = -1/6*v*(3^(1/2)*sin(3^(1/2)*k^(1/2)*t)-3*sin(k^(1/2)*t))/k^(1/2)

x2(t) = 1/6*v*(3^(1/2)*sin(3^(1/2)*k^(1/2)*t)+3*sin(k^(1/2)*t))/k^(1/2)

> simplify(subs(parms[2],eq15)); simplify(subs(parms[2],eq16));

x1(t) = -1/6*v*(3^(1/2)*sin(3^(1/2)*k^(1/2)*t)-3*sin(k^(1/2)*t))/k^(1/2)

x2(t) = 1/6*v*(3^(1/2)*sin(3^(1/2)*k^(1/2)*t)+3*sin(k^(1/2)*t))/k^(1/2)

基準モード

方程式 eq15およびeq16の初期条件がa1=0 あるいは a2=0 として指定されると面白い現象が起きます。

> p3:=unapply(subs([a1=0,a2=1,k=1,phi=0],rhs(eq15)),t);
p4:=unapply(subs([a1=0,a2=1,k=1,phi=0],rhs(eq16)),t);

p3 := proc (t) options operator, arrow; cos(t) end proc

p4 := proc (t) options operator, arrow; cos(t) end proc

> p5:=unapply(subs([a1=1,a2=0,k=1,phi=0],rhs(eq15)),t);
p6:=unapply(subs([a1=1,a2=0,k=1,phi=0],rhs(eq16)),t);

p5 := proc (t) options operator, arrow; cos(3^(1/2)*t) end proc

p6 := proc (t) options operator, arrow; -cos(3^(1/2)*t) end proc

基準モードのアニメーション

> oscillator(p3,p4,10,2*Pi/10);

[Plot]

> oscillator(p5,p6,10,2*Pi/sqrt(3)/10);

[Plot]

正規座標を使ってシステムを解きます。

変数の以下の様に変更してみましょう。

> eq19:=eta1(t)=x2(t)+x1(t); eq20:=eta2(t)=x2(t)-x1(t);

eq19 := eta1(t) = x2(t)+x1(t)

eq20 := eta2(t) = x2(t)-x1(t)

逆変換は次の様になります。

> solve({eq19,eq20},{x1(t),x2(t)});

{x1(t) = -1/2*eta2(t)+1/2*eta1(t), x2(t) = 1/2*eta1(t)+1/2*eta2(t)}

> eq21:=%[1];
eq22:=%%[2];

eq21 := x1(t) = -1/2*eta2(t)+1/2*eta1(t)

eq22 := x2(t) = 1/2*eta1(t)+1/2*eta2(t)

> eq3;
eq4;

(diff(x1(t), `$`(t, 2)))+2*k*x1(t)-k*x2(t) = 0

(diff(x2(t), `$`(t, 2)))-k*x1(t)+2*k*x2(t) = 0

正規座標に変換します。

> eq23:=simplify(subs([eq21,eq22],eq3)); eq24:=simplify(subs([eq21,eq22],eq4));

eq23 := (diff(-1/2*eta2(t)+1/2*eta1(t), `$`(t, 2)))-3/2*k*eta2(t)+1/2*k*eta1(t) = 0

eq24 := (diff(1/2*eta1(t)+1/2*eta2(t), `$`(t, 2)))+3/2*k*eta2(t)+1/2*k*eta1(t) = 0

上に示された2つの式を足し算および引き算します。結合された微分方程式が分離されたものになります。

> eq25:=eq24+eq23; eq26:=eq24-eq23;

eq25 := (diff(eta1(t), `$`(t, 2)))+k*eta1(t) = 0

eq26 := (diff(eta2(t), `$`(t, 2)))+3*k*eta2(t) = 0

分離された方程式を解きます。

> eq27:=dsolve({eq25,eta1(0)=0,D(eta1)(0)=v},{eta1(t)});

eq27 := eta1(t) = v*sin(k^(1/2)*t)/k^(1/2)

> eq28:=dsolve({eq26,eta2(0)=0,D(eta2)(0)=v},{eta2(t)});

eq28 := eta2(t) = 1/3*v*3^(1/2)*sin(3^(1/2)*k^(1/2)*t)/k^(1/2)

逆変換を行うと、前と同じ答えが得られます。

> simplify(subs([eq27,eq28],eq21));
simplify(subs([eq27,eq28],eq22));

x1(t) = 1/6*v*(-3^(1/2)*sin(3^(1/2)*k^(1/2)*t)+3*sin(k^(1/2)*t))/k^(1/2)

x2(t) = 1/6*v*(3*sin(k^(1/2)*t)+3^(1/2)*sin(3^(1/2)*k^(1/2)*t))/k^(1/2)