バネ付き振り子 - 武田のつぶや日記

色々悩んだけれどやっとバネ付き振り子っぽいものができた気がする。

方程式はソースを見ればわかるけど次のような感じ。
$\frac{dx}{dt} = v$

$\frac{d\theta}{dt} = \omega$

$\frac{dv}{dt} = (L + x)\omega^2 + gcos\theta - k\frac{x}{m}$

$\frac{d\omega}{dt} = -(2v\omega + gsin\theta)(L + x)$

ただし、xはバネの伸び、vは重りの速度、 $\theta$ は垂線とバネの角度、 $\omega$ は角速度、Lは自然長、gは重力加速度、kはバネ定数、mは重りの質量。

そしてソースはこんな感じ。

package lib;

import java.awt.Color;
import java.awt.Graphics;
import java.awt.image.BufferedImage;

import javax.swing.ImageIcon;
import javax.swing.JFrame;
import javax.swing.JLabel;

public class SamplePractice {

	static double gravity = 9.8;
	static double nuturalLength = 0.2;
	static double constantPendulum = 10;
	static double mass = 0.2;

	private static double funcIncrease(double t, double increase, double theata, double velocity, double omega){
		return velocity;
	}
	
	private static double funcTheta(double t, double increase, double theta, double velocity, double omega){
		return omega;
	}
	private static double funcVelocity(double t, double increase, double theta, double velocity, double omega){
		return (nuturalLength + increase) * omega * omega + gravity * Math.sin(theta) - constantPendulum * increase / mass;
	}
	private static double funcOmega(double t, double increase, double theta, double velocity, double omega){
		return -1 * (2 * velocity * omega + gravity * Math.sin(theta)) * (nuturalLength + increase);
	}
	
	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO 自動生成されたメソッド・スタブ

		BufferedImage bi = new BufferedImage(400, 300, BufferedImage.TYPE_INT_RGB);
		final Graphics g = bi.getGraphics();
		
		// 画面に表示
		JFrame frame = new JFrame("おもり付きバネ振り子");
		frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
		frame.setSize(400, 300);
		final JLabel label = new JLabel(new ImageIcon(bi));
		frame.add(label);
		frame.setVisible(true);
		
		// ここから
		new Thread(){
			public void run(){
				for(;;){
					try{
						Thread.sleep(300);
					}catch(InterruptedException e){
						// エラー処理
					}
					
					// 初期条件
					double h = 0.01; // 刻み幅
					double increase = 0;
					double theta = 0.5; // 弧度法(Math.PIを用いても、直接数値を入力しても問題なし。)
					double velocity = 0; // 初速度
					double omega = 0; // 初期角速度＝０＝ゆっくり手を離す
					
					int magnification = 400; // 拡大率(nuturalLengthを変更した際、球の様子が小さくて見えなくなる、あるいは大きすぎて見えない場合変更)					
					
					for(double t=0;t<40;t+=h){
						
						g.setColor(Color.WHITE);
						g.fillRect(0, 0, 400, 300);
						g.setColor(Color.BLACK);
						g.drawLine(200, 100, 203 + (int)(magnification*((nuturalLength + increase) * Math.sin(theta))), 103 + (int)(magnification*(nuturalLength * Math.cos(theta))));
						g.setColor(Color.RED);
						g.fillOval(200 + (int)(magnification*((nuturalLength + increase) * Math.sin(theta))), 100 + (int)(magnification*(nuturalLength * Math.cos(theta))), 7, 7);
						
						double increase1 = h * funcIncrease(t, increase, theta, velocity, omega);
						double theta1 = h * funcTheta(t, increase, theta, velocity, omega);
						double velocity1 = h * funcVelocity(t, increase, theta, velocity, omega);
						double omega1 = h * funcOmega(t, increase, theta, velocity, omega);
						
						double increase2 = h * funcIncrease(t + h/2, increase + increase1/2, theta + theta1/2, velocity + velocity1/2, omega + omega1/2);
						double theta2 = h * funcTheta(t + h/2, increase + increase1/2, theta + theta1/2, velocity + velocity1/2, omega + omega1/2);
						double velocity2 = h * funcVelocity(t + h/2, increase + increase1/2, theta + theta1/2, velocity + velocity1/2, omega + omega1/2);
						double omega2 = h * funcOmega(t + h/2, increase + increase1/2, theta + theta1/2, velocity + velocity1/2, omega + omega1/2);
						
						double increase3 = h * funcIncrease(t + h/2, increase + increase2/2, theta + theta2/2, velocity + velocity2/2, omega + omega2/2);
						double theta3 = h * funcTheta(t + h/2, increase + increase2/2, theta + theta2/2, velocity + velocity2/2, omega + omega2/2);
						double velocity3 = h * funcVelocity(t + h/2, increase + increase2/2, theta + theta2/2, velocity + velocity2/2, omega + omega2/2);
						double omega3 = h * funcOmega(t + h/2, increase + increase2/2, theta + theta2/2, velocity + velocity2/2, omega + omega2/2);
						
						double increase4 = h * funcIncrease(t + h, increase + increase3, theta + theta3, velocity + velocity3, omega + omega3);
						double theta4 = h * funcIncrease(t + h, increase + increase3, theta + theta3, velocity + velocity3, omega + omega3);
						double velocity4 = h * funcVelocity(t + h, increase + increase3, theta + theta3, velocity + velocity3, omega + omega3);
						double omega4 = h * funcOmega(t + h, increase + increase3, theta + theta3, velocity + velocity3, omega + omega3);
						
						double finalIncrease = (increase1 + 2*increase2 + 2*increase3 + increase4) / 6;
						double finalTheta = (theta1 + 2*theta2 + 2*theta3 + theta4) /6;
						double finalVelocity = (velocity1 + 2*velocity2 + 2*velocity3 + velocity4) / 6;
						double finalOmega = (omega1 + 2*omega2 + 2*omega3 + omega4) / 6;
						
						increase += finalIncrease;
						theta += finalTheta;
						velocity += finalVelocity;
						omega += finalOmega;
						
						label.repaint();
						
						try{
							Thread.sleep(10);
						}catch(InterruptedException e){
							// エラー処理
						}
					}
				}
			}
		}.start();

	}
}

あほらしい事でつまってたわ。
前回の単振動で、重力加速度を980[m/s^2]、自然長を100[m]として整合性？を保っていたのだけれど、よく考えたらこの方法だと今回のような式の場合(つまり、前回のように式がg/lみたいな感じになってなくてあるパラメータを何倍かした時に打ち消し合わないようなもの)うまくいくはずもなく、それでなぜ動かないんだとなっていた。
もっと言うと、本当は自然長0.5[m]、重力加速度9.8[m/s^2]とかでやりたいのだけれど、これをそのままやると自然長が短すぎてうまく動かなくなる。

なので今回は拡大率(magnification)という変数を用意して、実際に重りと糸を描画する際に得たパラメータを何百倍かして(ソースコードであれば400倍)、変化を大きく見せている。
最初は各初期値自体を何百倍かしていけばいいと単純に思っていたけれど、積み上げたときに誤差が出るので？こちらの方法で。

で、実際に動かしてみたものがこちら。

バネ定数を10[N/m]、重りの質量を0.2[Kg]、自然長は0.2[m]としたもの。
もう少しぬるぬる動けばいいなぁとか思ったけど、こんなものなのだろうか。
動きだけみても本当にあってるかどうか微妙なところだけれど、こんなものだろう。

ちなみにバネ定数を3[N/m]、重りの質量を0.5[Kg]、自然長は0.2[m]に値を変更したものはこちら。

荒ぶった動きをしています。

空気抵抗のない空間だとこんな感じの動きだろうか。
つまり、強いバネに軽い重りを取り付けると安定した動きをするようになり、弱いバネに重い重りを取り付けると不安定な動きをするようになると。

それからこの描画自体をどこでしているのかいまいちわかってなかったのだけれど、

BufferedImage で描画スペースを作り,
BufferedImageクラスのオブジェクトbiから描画するためのGraphicsクラスのオブジェクトgを作る。
新しいスレッドでgを用いて描画をして、
このBufferedImageクラスのオブジェクトbiをイメージアイコン→JLabelとして取り込み
そのラベルをframeにaddしている。

という感じなのかな。
知識の浅い自分には新鮮な書き方だった。
なんとなく。

でもやり方思いついて本当に良かった。
考えながら寝たら意外と朝起きたらこうすればいいんじゃね？って思うものだね。