for構文で算出した​複数の解を一つの行列​の行に羅列するにはど​うすればよいでしょう​か。

9 visualizaciones (últimos 30 días)
m17td024
m17td024 el 21 de Nov. de 2018
Comentada: m17td024 el 28 de Nov. de 2018
for構文において算出した複数(非常に多い)の解を一つに行列の行に羅列するには、どうすればよいのでしょうか。
syms Z1 Z2;
for i=0:0.01:3.9
Z1=1+i;
a1=(Z1*sin(Z1)-Z1^2*cos(Z1))/(2-2*cos(Z1)-Z1*sin(Z1));
b1=(Z1^2-Z1*sin(Z1))/(2-2*cos(Z1)-Z1*sin(Z1));
a2=(Z2*sin(Z2)-Z2^2*cos(Z2))/(2-2*cos(Z2)-Z2*sin(Z2));
b2=(Z2^2-Z2*sin(Z2))/(2-2*cos(Z2)-Z2*sin(Z2)); %安定関数
K1=2*((a2-b2)+3/2)+(3/4+12);
K2=2*((a1-b1)+3/2)+(3/4+12);
K3=2*((a1-b1+a2-b2)+(3/4+4));
K5=-(3/4+2); %行列要素
K=[K3 0 0 0 0 0;0 K3 0 0 0 0;0 0 K1 K5 0 0;0 0 K5 K2 0 0;...
0 0 0 0 K1 K5;0 0 0 0 K5 K2]; %座屈モードⅠの行列
ka=K(1,1); %モードタイプaの行列
solx=vpasolve(det(ka)==0,Z2,[0 10]);
end
solxに算出した解が随時上書きされてしまうことも問題点です。
一つの行列Aに、
A=[i=0のときの解;i=0.01のときの解;・・・;i=3.9のときの解]
のようにしたいのですが。
よろしくお願い致します。

Respuesta aceptada

Yoshio
Yoshio el 21 de Nov. de 2018
Editada: Yoshio el 21 de Nov. de 2018
一つの行列が、一列のベクトルという意味でしたら、以下のコードではいかがでしょうか?
syms Z1 Z2;
n = 0;
for i=0:0.01:3.9
n = n+1;
Z1=1+i;
a1=(Z1*sin(Z1)-Z1^2*cos(Z1))/(2-2*cos(Z1)-Z1*sin(Z1));
b1=(Z1^2-Z1*sin(Z1))/(2-2*cos(Z1)-Z1*sin(Z1));
a2=(Z2*sin(Z2)-Z2^2*cos(Z2))/(2-2*cos(Z2)-Z2*sin(Z2));
b2=(Z2^2-Z2*sin(Z2))/(2-2*cos(Z2)-Z2*sin(Z2)); %安定関数
K1=2*((a2-b2)+3/2)+(3/4+12);
K2=2*((a1-b1)+3/2)+(3/4+12);
K3=2*((a1-b1+a2-b2)+(3/4+4));
K5=-(3/4+2); %行列要素
K=[K3 0 0 0 0 0;0 K3 0 0 0 0;0 0 K1 K5 0 0;0 0 K5 K2 0 0;...
0 0 0 0 K1 K5;0 0 0 0 K5 K2]; %座屈モードⅠの行列
ka=K(1,1); %モードタイプaの行列
solx=vpasolve(det(ka)==0,Z2,[0 10]);
a(n) = solx; %行ベクトルに保存
end
A = a.' % 解が複素数の場合を考慮
  1 comentario
m17td024
m17td024 el 22 de Nov. de 2018
列ベクトルとして、解を羅列することができました。
誠にありがとうございます。

Iniciar sesión para comentar.

Más respuestas (3)

Yoshio
Yoshio el 24 de Nov. de 2018
シンボリックでも指定の大きさでメモリの事前割り当てができますので、修正してみました。
syms Z1 Z2;
index = 0:0.01:3.9;
m = length(index);
A = sym(zeros(m,1));
n = 0;
for i = 0:0.01:3.9;
n = n+1;
Z1=1+i;
a1=(Z1*sin(Z1)-Z1^2*cos(Z1))/(2-2*cos(Z1)-Z1*sin(Z1));
b1=(Z1^2-Z1*sin(Z1))/(2-2*cos(Z1)-Z1*sin(Z1));
a2=(Z2*sin(Z2)-Z2^2*cos(Z2))/(2-2*cos(Z2)-Z2*sin(Z2));
b2=(Z2^2-Z2*sin(Z2))/(2-2*cos(Z2)-Z2*sin(Z2)); %安定関数
K1=2*((a2-b2)+3/2)+(3/4+12);
K2=2*((a1-b1)+3/2)+(3/4+12);
K3=2*((a1-b1+a2-b2)+(3/4+4));
K5=-(3/4+2); %行列要素
K=[K3 0 0 0 0 0;0 K3 0 0 0 0;0 0 K1 K5 0 0;0 0 K5 K2 0 0;...
0 0 0 0 K1 K5;0 0 0 0 K5 K2]; %座屈モードⅠの行列
ka=K(1,1); %モードタイプaの行列
%solx=vpasolve(det(ka)==0,Z2,[0 10]);
A(n) = vpasolve(det(ka)==0,Z2,[0 10]);
end
  1 comentario
m17td024
m17td024 el 26 de Nov. de 2018
ご改善ありがとうございます。
もう一つお聞きしてもよろしいでしょうか。
Z2の計算において、0から10の範囲で解が複数個存在するのですが、
その中から最小解を算出するにはどうすればよいでしょうか。

Iniciar sesión para comentar.


Yoshio
Yoshio el 26 de Nov. de 2018
Editada: Yoshio el 26 de Nov. de 2018
解の精度をどの程度求めるかにもよると思いますが、シンボリック関数 ka を
ka=K(1,1); %モードタイプaの行列
y = matlabFunction(det(ka));
とすることで、通常の関数になりますから、一変数の最小化問題として、fzero等を用いて
[x,fval]= fzero(y,x0)
のように解とその時の評価値を調べ、xが0から10の範囲でabs(fval)が最小となるxを求めれば良いと思います。非線形の問題ですので、初期値x0と解x に関する許容誤差 の設定は、適切に行う必要があります。
  1 comentario
m17td024
m17td024 el 28 de Nov. de 2018
お教えくださったもので計算すると、単体で解を出すと最小値を得ることができましたが、
ループ構文の中に当てはめてみると、途中から解がすべて0になってしまいます。
具体的には、Z1が1からスタートして、初めの10回目のループまでは、順調に最小解が算出されますが、そのあとのループにおいては、すべて0となってしまいます。
単体でZ1=3の場合で計算すると、解がででくるのですが…
問題点があれば、お教えください。
syms Z1 Z2;
index = 0:0.01:3.9;
m = length(index);
A = sym(zeros(m,1));
n = 0;
for i = 0:0.1:3.9;
n = n+1;
Z1=1+i;
a1=(Z1*sin(Z1)-Z1^2*cos(Z1))/(2-2*cos(Z1)-Z1*sin(Z1));
b1=(Z1^2-Z1*sin(Z1))/(2-2*cos(Z1)-Z1*sin(Z1));
a2=(Z2*sin(Z2)-Z2^2*cos(Z2))/(2-2*cos(Z2)-Z2*sin(Z2));
b2=(Z2^2-Z2*sin(Z2))/(2-2*cos(Z2)-Z2*sin(Z2)); %安定関数
K1=2*((a2-b2)+3/2)+(3/4+12);
K2=2*((a1-b1)+3/2)+(3/4+12);
K3=2*((a1-b1+a2-b2)+(3/4+4));
K5=-(3/4+2); %行列要素
K=[K3 0 0 0 0 0;0 K3 0 0 0 0;0 0 K1 K5 0 0;0 0 K5 K2 0 0;...
0 0 0 0 K1 K5;0 0 0 0 K5 K2]; %座屈モードⅠの行列
ka=K(1,1); %モードタイプaの行列
y = matlabFunction(det(ka));
x0=1;
A(n)= fzero(y,x0);
end

Iniciar sesión para comentar.


Yoshio
Yoshio el 28 de Nov. de 2018
Editada: Yoshio el 28 de Nov. de 2018
まずyをxの関数としてプロットし、初期値の選択が適切か検討することから始めましょう。一つの固定したi について、x0を[0 10]の範囲で振ってfzero(y,x0);の解と、精度を求めてそれから適切な解を選択して、A(n)に入れてはどうでしょうか。
  1 comentario
m17td024
m17td024 el 28 de Nov. de 2018
途中式でミスがあり、修正して関数をプロットすると、解は複数個存在しませんでした。
そのため、最初にご回答いただいたもので問題が解決しそうです。
お手数お掛けし、誠に申し訳ございませんでした。
ご親切にご回答いただきありがとうございました。

Iniciar sesión para comentar.

Categorías

Más información sobre 数学 en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!