plotができない

30 views (last 30 days)
kazuma kaneda
kazuma kaneda on 21 Nov 2022
Moved: Atsushi Ueno on 25 Nov 2022
微分方程式をルンゲクッタ法を用いて解こうとしています.
添付ファイルのコードを実行しても.グラフを表示できませんでした.
for loopの部分がおかしいと思うのですが,それがどこなのかわかりません.

Answers (2)

Atsushi Ueno
Atsushi Ueno on 21 Nov 2022
> 添付ファイルのコードを実行しても.グラフを表示できませんでした
グラフは表示されています。分かり易い様、間に線も引いてみました。
問題はプログラムを実行し終えてプロットした 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')

kazuma kaneda
kazuma kaneda on 24 Nov 2022
Moved: Atsushi Ueno on 25 Nov 2022
回答ありがとうございます.
x についてですが,配列ではなくスカラーで表したいです.
後半は解決しました.
  1 Comment
Atsushi Ueno
Atsushi Ueno on 24 Nov 2022
Moved: Atsushi Ueno on 25 Nov 2022
> 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

Sign in to comment.

Categories

Find more on ループと条件付きステートメント in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!