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

1 次查看（过去 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)
のように解とその時の評価値を調べ、ｘが0から10の範囲でabs(fval)が最小となるｘを求めれば良いと思います。非線形の問題ですので、初期値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 Center 和 File Exchange 中查找有关 数学 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!