plotができない

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

답변(2개)

Atsushi Ueno
Atsushi Ueno 2022년 11월 21일
> 添付ファイルのコードを実行しても.グラフを表示できませんでした
グラフは表示されています。分かり易い様、間に線も引いてみました。
問題はプログラムを実行し終えてプロットした 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 2022년 11월 24일
Moved: Atsushi Ueno 2022년 11월 25일
回答ありがとうございます.
x についてですが,配列ではなくスカラーで表したいです.
後半は解決しました.
  댓글 수: 1
Atsushi Ueno
Atsushi Ueno 2022년 11월 24일
Moved: Atsushi Ueno 2022년 11월 25일
> 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

댓글을 달려면 로그인하십시오.

범주

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

제품


릴리스

R2021b

Community Treasure Hunt

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

Start Hunting!