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

v=A*omega*cos(omega*t+alpha);

上記の単振動の公式を使ってこれをオイラー法で近似したいと思いますが、

for ループがうまくいかないです。

恐らくｖをスカラー値にしたらできると思うのですが、これではｙがスカラー値になってしまっていると思うのですが

どのように変更するのかわかりません。

もし分かれば教えていただきたいです

よろしくお願いいたします。

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');

Answer 1

[Takumi]

エラーが出る理由は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');