単振動をオイラー法で数値シミュレーション

63 views (last 30 days)
sz
sz on 15 May 2020
Answered: Takumi on 16 May 2020
v=A*omega*cos(omega*t+alpha);
上記の単振動の公式を使ってこれをオイラー法で近似したいと思いますが、
for ループがうまくいかないです。
恐らくvをスカラー値にしたらできると思うのですが、これではyがスカラー値になってしまっていると思うのですが
どのように変更するのかわかりません。
もし分かれば教えていただきたいです
よろしくお願いいたします。
A=10;
alpha=pi./2;
omega=2*pi;
h=1/250;
y(1)=10;
N=1+4/h;
t=0:h:4;
for i=1:1:N;
v=A*omega*cos(omega*t+alpha);
v0=0;
y(i+1)=y(i)+v*h;  %ここでエラーが出ます。
t(i+1)=t(i)+h;
end
yy=y(1,1:N);
plot(t(1,:),y(1,:),'-o');
hold on
plot(t,v,'-o');

Accepted Answer

Takumi
Takumi on 16 May 2020
エラーが出る理由はy(i+1)という一つの要素にy(i)+v*hという配列を代入しようとしているからです。y(i)+v*hが配列になってしまうのはvが配列であるからで、vが配列になってしまうのは
v=A*omega*cos(omega*t+alpha);
とすべてのtに対して計算しているからです。やりたいことはおそらく
v=A*omega*cos(omega*t(i)+alpha);
ではないでしょうか。そして他にも初歩的なミスがあるので例えば以下のようにするといいと思います。参考にしてください。
A=10; % 振幅
alpha=pi./2; % 位相
omega=2*pi; % 角振動数
h=1/250; % 時間刻み
N=1+4/h; % 計算点数
% 初期条件
t(1)=0; % 初期時間
y(1)=10; % 初期位置
v(1)=A*omega*cos(omega*t(1)+alpha); % 初速度
for i=1:N-1
t(i+1)=t(i)+h; % 時間
v(i+1)=A*omega*cos(omega*t(i+1)+alpha); % 速度
y(i+1)=y(i)+v(i)*h; % 位置
end
plot(t,y,'-o');
hold on
plot(t,v,'-o');
また、自分だったらこのように書きます。
A = 10; % 振幅
alpha = pi./2; % 位相
omega = 2*pi; % 角振動数
h = 1/250; % 時間刻み
t = 0:h:4; % 時間配列
v = A*omega*cos(omega*t+alpha); % 速度
% 初期条件
y = zeros(size(t)); % 配列初期化
y(1) = 10; % 初期位置
for i = 1:length(t)-1
y(i+1) = y(i)+v(i)*h; % 位置
end
plot(t,y,'-o');
hold on
plot(t,v,'-o');

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!