plotができない

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

Respuestas (2)

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

Iniciar sesión para comentar.

Categorías

Más información sobre ループと条件付きステートメント en Help Center y File Exchange.

Productos


Versión

R2021b

Community Treasure Hunt

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

Start Hunting!