基本的な力学問題2

3. 基本的な力学問題 (続き)
2014/10/29
1
1.1
非線形振動子
運動方程式
前回の実験で、調和振動子の数値シミュレーションを行いました。調和振動子
とは、微分方程式
d2 x
(1.1)
m 2 = −kx
dt
で表される振動子のことです。
振幅が大きくなると、右辺が振幅 x に比例するだけでなく、振幅の高次の項に
も依存することが予想されます。そのような、高次の項を「非線形項」と呼び、非
線形項を含む振動子を「非線形振動子」と呼びます。非線形振動子の運動方程式
として、二つの定数 k と λ を用いて
m
d2 x
= −kx − λx3
2
dt
(1.2)
を考えることにします。この運動方程式を二つの変数 (x, v) の連立微分方程式に分
けると以下のようになります。
m
dv
= −kx − λx3
dt
dx
= v
dt
(1.3)
(1.4)
前回配布の Java ソースコードでは、調和振動子のクラス HarmonicOscillator
の場合だけに、その運動方程式が実装されていました。非線形振動子のクラス
NonlinearOscillator では、そのコンストラクタが以下のように運動方程式部
分が実装されていません (Program 1.1)。ここで、引数 k は k/m に相当し、lambda
は λ/m に相当します。ここに、式 (1.2) に相当する部分を記述しましょう。
1
public NonlinearOscillator(double x, double v,
double k, double lambda) {
super(2, x, v);
this.k = k;
this.lambda = lambda;
equation = (double xx, double[] yy) -> {
double dy[] = new double[numVar];
return dy;
};
}
Program 1.1: NonlinearOscillator のコンストラクタ部
1.2
実行
非線形振動子のクラス NonlinearOscillator の main メソッドは、Program 1.2
となっています。出力ファイルは output.txt です。前回のファイルと区別できる
ように、nonlinear.txt に変更し、実行しましょう。ファイル nonlinear.txt が
出来ていることを確認します。
出力ファイルに対応して作図する gnuplot 用のスクリプトを Program 1.3 に示
します。ここで a と b は、振動の上端と下端を示すための値を適当に入れていま
す。シミュレーションで使用したパラメタ m = 1、k = 1、λ = 4 に対応した値を
設定して、振動の上端と下端が正しい事を確かめましょう。
2
public static void main(String[] args) throws IOException {
NonlinearOscillator sys =
new NonlinearOscillator(0., 1., 1., 4.);
double t = 30.;
int nstep = 100000;
List<Point2D.Double> points = sys.evolution(t, nstep);
BufferedWriter out = FileIO.openWriter("output.txt");
for (Point2D.Double p : points) {
FileIO.writeSSV(out, p.x,p.y);
}
}
Program 1.2: NonlinearOscillator の main 部
set terminal png enhanced
set xlabel "t"
set ylabel "x"
set yrange [-1.5:1.5]
set xrange [0:20]
set ytic 0.5
set title "Nonlinear Oscillator"
set output "NonlinearOscillator.png"
a = 1
b = 1
plot "output.txt", a lw 3, b lw 3
Program 1.3: nonlinear.plt
3
public DumpedOscillator(double x, double v,
double k, double lambda) {
super(2, x, v);
this.k = k;
this.lambda = lambda;
equation = (double xx, double[] yy) -> {
double dy[] = new double[numVar];
return dy;
};
}
Program 2.1: DumpedOscillator のコンストラクタ
2
減衰振動
運動方程式の右辺に速度に比例した項を入れることで摩擦を表すことができます。
m
d2 x
= −kx − γv
dt2
(2.1)
前回配布の DumpedOscillator のコンストラクタを完成させ、実行して、振幅が
時間とともに小さくなることを示しましょう。
4

Download Report