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

6 visualizzazioni (ultimi 30 giorni)
m17td024
m17td024 il 21 Nov 2018
Commentato: m17td024 il 28 Nov 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のときの解]
のようにしたいのですが。
よろしくお願い致します。

Risposta accettata

Yoshio
Yoshio il 21 Nov 2018
Modificato: Yoshio il 21 Nov 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 Commento
m17td024
m17td024 il 22 Nov 2018
列ベクトルとして、解を羅列することができました。
誠にありがとうございます。

Accedi per commentare.

Più risposte (3)

Yoshio
Yoshio il 24 Nov 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 Commento
m17td024
m17td024 il 26 Nov 2018
ご改善ありがとうございます。
もう一つお聞きしてもよろしいでしょうか。
Z2の計算において、0から10の範囲で解が複数個存在するのですが、
その中から最小解を算出するにはどうすればよいでしょうか。

Accedi per commentare.


Yoshio
Yoshio il 26 Nov 2018
Modificato: Yoshio il 26 Nov 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 Commento
m17td024
m17td024 il 28 Nov 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

Accedi per commentare.


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

Accedi per commentare.

Community Treasure Hunt

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

Start Hunting!