matlab 衝突のシミュレーション

質問

現在、Matlabで簡単な衝突のシミュレーションを行っていますが、うまくいかず、どなたがご教授ください。

【やりたいこと】
ボールを自由落下させた時の地面との衝突のシミュレーションを、
Simulinkではなく、自分で作成したmfileから実現させたい。

【質問】
以下のプログラムでodeの計算が終わらず、どのようにすればいいのか
ご教授ください。

=======================
【diff_freefall.m】
function dx = diff_freefall(t,x0)
g = 9.81; % 重力
k = 1; % 反発係数

A = [0 1;
0 0];

B = [ 0;
-g];

% 衝突判定
if x0(1)<=0 && x0(2)<0
x0(1) = 0;
x0(2) = -k*x0(2);
end

dx = A * x0 + B;
=======================

=======================
【simu.m】
Ts = 0; %シミュレーション開始時間
Te = 10; %シミュレーション終了時間
y0 = [10; 15;]; %初期位置、初期速度

[T,Y] = ode45('diff_freefall',[Ts:0.1:Te],y0);

subplot(2,1,1);
plot(T,Y(:,1),'g');
subplot(2,1,2);
plot(T,Y(:,2),'r');
=======================

【補足】
・以下のmathworksのサイトのサンプル(simulink版)を、mfileで記述したいです。
<http://jp.mathworks.com/help/simulink/examples/simulation-of-a-bouncing-ball_ja_JP.html>

・すべてのodeのタイプを試しましたが、ずっと終わらない、
　もしくは途中でシミュレーションが止まります。

・衝突した時に、odeの積分の累積値がクリアできていないことが原因だと思うのですが、
　クリアする方法がわからず困っています。

回答No.2

jjorz · ベストアンサー

Scilabでやってみました。
Matlabに書き直すのはそんなに難しくないはずです。

clear;

// *** 計算条件 ***
g = 9.81; // 重力加速度
k = 0.8; // 反発係数
// 時間
ts = 0; // 開始時刻
te = 10; // 終了時刻
// 初期条件
x0 = 10; // 初期位置
v0 = 15; // 初速度

// *** 常微分方程式の定義 ***
function dx = fall(t,x)
dx(1) = x(2); // dx/dt = v
dx(2) = -g; // dv/dt = -g
endfunction

// *** 解くべき非線形方程式 ***
function x = bound(t)
X = ode([x0; v0], ts, t, fall);
x = X(1);
endfunction

// *** メイン ***
// 初期化
tb = ts; // 衝突時刻
X = []; // プロット用の位置と速度
T = []; // プロット用の時間
// 終了時刻まで繰り返し
while tb < te
// 衝突時刻を計算
tb = fsolve(2 * te, bound);
// プロット用の時間ベクトルを作成
if tb > te then
// 終了時刻まで
time = linspace(ts, te);
else
// 衝突時刻まで
time = linspace(ts, tb);
end
T = [T, time];
// プロット用の位置と速度を計算
X = [X, ode([x0; v0], ts, time, fall)];
// 次の繰返しのための初期条件
x0 = 0; // 初期位置
v0 = -k * X(2,$); // 初速度
ts = tb; // 計算開始時刻
end

// *** グラフのプロット ***
// 位置のプロット
subplot(2,1,1);
plot(T, X(1,:),'g');
xgrid(color("gray"));
xlabel("Time");
ylabel("Position");
// 速度のプロット
subplot(2,1,2);
plot(T, X(2,:),'r');
xgrid(color("gray"));
xlabel("Time");
ylabel("Velocity");

お礼

理解することができました！
コードまで作成していただき、ありがとうございます。

回答No.1

jjorz

衝突ごとに別々にodeを計算してはどうですか？

1. 微分方程式は衝突判定を含めない
2. odeとfsolveを組み合わせて衝突までの計算をする
3. 反発後の速度を初速度として、次の衝突までode+fsolve
4. 2と3を終了時刻まで繰り返し

お礼

回答ありがとうございます。

==========================================
>> 衝突ごとに別々にodeを計算してはどうですか？
→　やりたいことはまさにこれです。

>>1. 微分方程式は衝突判定を含めない
　→　理解しました。

>> 2. odeとfsolveを組み合わせて衝突までの計算をする
　→　"fsolve"、知りませんでした。

>> 3. 反発後の速度を初速度として、次の衝突までode+fsolve
　→　”次の衝突までode+fsolve”が理解できていません。

>> 4. 2と3を終了時刻まで繰り返し
　→　理解しました。
==========================================

これからfsolveを調べてみますが、odeで微分方程式を計算中、
時間以外の条件で終了させることができるということですよね？
引き続き、がんばります。