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

조회 수: 23 (최근 30일)
m17td024
m17td024 2018년 11월 21일
댓글: m17td024 2018년 11월 28일
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のときの解]
のようにしたいのですが。
よろしくお願い致します。

채택된 답변

Yoshio
Yoshio 2018년 11월 21일
편집: Yoshio 2018년 11월 21일
一つの行列が、一列のベクトルという意味でしたら、以下のコードではいかがでしょうか?
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
m17td024
m17td024 2018년 11월 22일
列ベクトルとして、解を羅列することができました。
誠にありがとうございます。

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

추가 답변 (3개)

Yoshio
Yoshio 2018년 11월 24일
シンボリックでも指定の大きさでメモリの事前割り当てができますので、修正してみました。
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
m17td024
m17td024 2018년 11월 26일
ご改善ありがとうございます。
もう一つお聞きしてもよろしいでしょうか。
Z2の計算において、0から10の範囲で解が複数個存在するのですが、
その中から最小解を算出するにはどうすればよいでしょうか。

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


Yoshio
Yoshio 2018년 11월 26일
편집: Yoshio 2018년 11월 26일
解の精度をどの程度求めるかにもよると思いますが、シンボリック関数 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
m17td024
m17td024 2018년 11월 28일
お教えくださったもので計算すると、単体で解を出すと最小値を得ることができましたが、
ループ構文の中に当てはめてみると、途中から解がすべて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

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


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

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

카테고리

Help CenterFile Exchange에서 数学에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!