解決済み

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

  • すぐに回答を!
  • 質問No.8885855
  • 閲覧数2030
  • ありがとう数4
  • 気になる数2
  • 回答数2
  • コメント数0

お礼率 100% (2/2)

現在、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

ベストアンサー率 66% (2/3)

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");
お礼コメント
dankichi1218

お礼率 100% (2/2)

理解することができました!
コードまで作成していただき、ありがとうございます。
投稿日時 - 2015-01-14 21:30:53

その他の回答 (全1件)

  • 回答No.1

ベストアンサー率 66% (2/3)

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

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

お礼率 100% (2/2)

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

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

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

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

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

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

これからfsolveを調べてみますが、odeで微分方程式を計算中、
時間以外の条件で終了させることができるということですよね?
引き続き、がんばります。
投稿日時 - 2015-01-11 12:33:08
結果を報告する
このQ&Aにはまだコメントがありません。
あなたの思ったこと、知っていることをここにコメントしてみましょう。
関連するQ&A
AIエージェント「あい」

こんにちは。AIエージェントの「あい」です。
あなたの悩みに、OKWAVE 3,600万件のQ&Aを分析して最適な回答をご提案します。

その他の関連するQ&A、テーマをキーワードで探す

キーワードでQ&A、テーマを検索する

ピックアップ

ページ先頭へ