plotができない

微分方程式をルンゲクッタ法を用いて解こうとしています．

添付ファイルのコードを実行しても．グラフを表示できませんでした．

for loopの部分がおかしいと思うのですが，それがどこなのかわかりません．

Answer 1

[Atsushi Ueno]

> 添付ファイルのコードを実行しても．グラフを表示できませんでした

グラフは表示されています。分かり易い様、間に線も引いてみました。

問題はプログラムを実行し終えてプロットした x = [ 0, 0.0559 ] の値が意図通りかどうかです。

> for loopの部分がおかしいと思うのですが，それがどこなのかわかりません

ロジックは確認していません。

for loopで10回ループを回しているのに対し、変数 iter の値がずっと1のままになっています。

そのまま x(iter+1) や g(iter+1) の値を計算しても、10回上書きするだけになります。

clc; clear all;
% 係数の決定
m=10; c=20; k=25; mu=5; e=0.25; w=10;
ep=c/m;
mr=m/mu;
wn=k/m;
% 積分区間の設定
t1=0.0;
t2=1.0;
% 初期値
x1=0.0; g1=0.0;
% きざみ
n=10;
%
delt=(t2-t1)/n;
%
iter=1;
x(iter)=x1;
t(iter)=t1;
g(iter)=g1;
for i=1:n
    k1x=delt.*g(iter);
    k1g=delt.*(-2*ep*g(iter)-wn.^2*x(iter)+mr*e*w.^2*sin(w*t(iter)));

    k2x=delt.*(g(iter)+k1g./2);
    k2g=delt.*(-2*ep*(g(iter)+k1g./2)-wn.^2*(x(iter)+k1x./2)+mr*e*w^2*sin(w*(t(iter)+delt/2)));

    k3x=delt.*(g(iter)+k2g./2);
    k3g=delt.*(-2*ep*(g(iter)+k2g./2)-wn.^2*(x(iter)+k2x./2)+mr*e*w.^2*sin(w*(t(iter)+delt/2)));

    k4x=delt.*(g(iter)+k3g/2);
    k4g=delt.*(-2*ep*(g(iter)+k3g)-wn.^2*(x(iter)+k3g)+mr*e*w.^2*sin(t(iter)+delt));

    x(iter+1)=x(iter)+(k1x+2.*k2x+2.*k3x+k4x)./6;
    g(iter+1)=g(iter)+(k1g+2.*k2g+2.*k3g+k4g)./6;
end
plot(x,'r-o')

Answer 2

[kazuma kaneda]

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

x についてですが，配列ではなくスカラーで表したいです．

後半は解決しました．

[Atsushi Ueno]

> x についてですが，配列ではなくスカラーで表したいです．

では plot 関数を for ループの内側に入れて都度描画していくのはどうでしょう。線で繋げる為、一度に前回値から今回値までを描画して重ねていきます。

ループ内で t(三角関数の時間と思われるパラメータ) も動かないとならないような気がしますが何を計算しているの良くわかりません。

m=10; c=20; k=25; mu=5; e=0.25; w=10; ep=c/m; mr=m/mu; wn=k/m;
t1=0.0; t2=1.0; x1=0.0; g1=0.0; n=10; delt=(t2-t1)/n;
x=x1;
t=t1;
g=g1;
figure; % プロット軸を収めるfigure画面を出す
hold on % プロット描画内容を保持する(再描画時にクリアしない)
for i=1:n
    k1x=delt.*g;
    k1g=delt.*(-2*ep*g-wn.^2*x+mr*e*w.^2*sin(w*t));
    k2x=delt.*(g+k1g./2);
    k2g=delt.*(-2*ep*(g+k1g./2)-wn.^2*(x+k1x./2)+mr*e*w^2*sin(w*(t+delt/2)));
    k3x=delt.*(g+k2g./2);
    k3g=delt.*(-2*ep*(g+k2g./2)-wn.^2*(x+k2x./2)+mr*e*w.^2*sin(w*(t+delt/2)));
    k4x=delt.*(g+k3g/2);
    k4g=delt.*(-2*ep*(g+k3g)-wn.^2*(x+k3g)+mr*e*w.^2*sin(t+delt));
    oldx = x; oldg = g;
    x=x+(k1x+2.*k2x+2.*k3x+k4x)./6;
    g=g+(k1g+2.*k2g+2.*k3g+k4g)./6;
    plot([i-1 i],[oldx x],'r-o',[i-1 i],[oldg g],'b-*');
    % tが動かないのでiをx軸に、xとgをy軸に設定したplotを表示
end