単振動をオイラー法で数値シミュレーション
显示 更早的评论
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');
采纳的回答
Takumi
2020-5-16
エラーが出る理由は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');
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 パフォーマンスとメモリ 的更多信息
