ある条件をクリアする​まで、演算を繰り返す​にはどうすればよいで​しょうか。

以下のプログラミングにおいて、条件をクリアするまで演算を繰り返すようにするには、どのようにすればよいのでしょうか。

変数Nを固定し、Qを15~で小数第三位で増分していき(15,15.001,15.002…)、Aがもっとも0に近似した段階のQを求めるという形にしたいのです。（方程式でQをダイレクトに求める方法(A=0の)では、式が複雑すぎて、途中でフリーズしてしまいます。）

以下のプログラミングでは、15~16でQを増分し、算出されたAのなかからもっとも0に近いものを手作業で見つけ、その時のQを採用するというやり方になり、非常に困っております。

どうぞよろしくお願い致します。

syms N Q; 
N=2*pi^2;
index = 0:0.001:1;
m = length(index);
A = sym(zeros(m,1));
n = 0;
for i = 0:0.001:1;

    n = n+1;

    Q=15+i;

Z1=(N+2^0.5*Q/2)^0.5;
Z2=(-(-N+2^0.5*Q/2))^0.5;
Z3=(N-2^0.5*Q/2)^0.5;
Z4=(-(-N-2^0.5*Q/2))^0.5;

ah1=(Z1*sinh(Z1)-Z1^2*cosh(Z1))/(2*cosh(Z1)-2-Z1*sinh(Z1));
bh1=(Z1^2-Z1*sinh(Z1))/(2*cosh(Z1)-2-Z1*sinh(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));
ah3=(Z3*sinh(Z3)-Z3^2*cosh(Z3))/(2*cosh(Z3)-2-Z3*sinh(Z3));
bh3=(Z3^2-Z3*sinh(Z3))/(2*cosh(Z3)-2-Z3*sinh(Z3));
a4=(Z4*sin(Z4)-Z4^2*cos(Z4))/(2-2*cos(Z4)-Z4*sin(Z4));
b4=(Z4^2-Z4*sin(Z4))/(2-2*cos(Z4)-Z4*sin(Z4)); 

K1=2*(2+3/2)+(3/4+3*(ah1+a2+ah3+a4)/4);
K2=2*(2+3/2)+(3/4+3*(ah1+a2+ah3+a4)/4);
K3=8+3/2+(ah1+a2+ah3+a4)/2;
K4=-(ah1-a2+ah3-a4)/4;
K5=(ah1-a2-ah3+a4)/(2*2^0.5);
K6=3*(bh1-b2+bh3-b4)/4;
K7=-3/4-(bh1+b2+bh3+b4)/4;
K8=(bh1+b2-bh3-b4)/(2*2^0.5);
K9=(ah1+a2-ah3-a4)/(2*2^0.5);
K10=(bh1-b2-bh3+b4)/(2*2^0.5);
K11=(bh1-b2+bh3-b4)/2; 

K=[K1 K4 K5 K6 K7 K8;K4 K2 K9 K7 K6 K10;K5 K9 K3 K8 K10 K11;...
    K6 K7 K8 K1 K4 K5;K7 K6 K10 K4 K2 K9;K8 K10 K11 K5 K9 K3]; 

A(n) =det(K);
end

Answer 1

[Kazuya]

条件をクリアするまで実施する方法、ではありませんが fminbnd関数がまさに目的にあっているような気がします。f(Q) = abs(A - 0) が最小値となる Q を 15から16の間から見つけることができるかもしれません。Q の値によってAの値が変わってくるのかによってうまくいかないこともありますが。

https://jp.mathworks.com/help/matlab/ref/fminbnd.html

サンプルコードがいくつかあるので参考にしてみてください。

[m17td024]

ご回答ありがとうございます。

今回の場合はQについての一変数関数であるため、教えていただいたfminbnd関数で、

理論上は算出可能なのですが、A=det(K)がQについての非常に複雑な関数であるため、

計算に非常に時間を要するという問題点(A=det(K)=0の方程式を解こうとすると、

少なくとも1時間では結果が出ない)があり、それに対応するために、実際にQに数値を地道に

代入して解きたいのですが…

実際にQに数値を代入して、結果を得る方法はございませんでしょうか。

[Kazuya]

Q に数値を代入するということは、 A = f(Q) : A が Q の関数として定義できて、この関数に対して 数値（Qの値）を代入することで A の値を計算したいということですか？

それぞれのQにおける A の値を計算して、A がお望みの値（最も小さい値）になる Q の値を見つけるというのが目的、、と理解していたのですが違いますでしょうか。。なんどもすいません。

[m17td024]

すいません。私の勘違いでした。おっしゃる通りです。

代入した後のAについて、fminbnd関数を適用させるということでよろしいのでしょうか。

その場合、Qの値が変わるごとにAが算出され、それを連続的にfminbnd関数に適用させるには、どうすればよいのでしょうか。

[Kazuya]

もし、、いろいろなQでのAの値も求めたいのであれば話は別ですが、Aが最も0に近くなる時のQだけを求めるなら fminbnd を使う方法が楽なのでは、という提案でした。。計算時間も短く済むと思います。

A = f(Q) となる関数をまず作ってください。作った関数を fun.m とすると

x = fminbnd(fun,15,16)

で fun が最小の値を出すときの Q の値を計算します。

「fun が最小」なので、「Aの値が最も0に近い」ものを求めたければ、fun は abs(A) と絶対値を返すような関数である必要がありますね。

言葉で説明するより、 fminbnd 関数のページの例題を眺めてイメージをつかんでもらう方が早い気がしますが、例えば

functoion absA = fun(Q)
(ここにAを計算するプログラム)
absA = abs(A);
end

てな感じになると思いますが、、どうでしょうか？

[m17td024]

"A = f(Q) となる関数"ですが、そもそもQの変数を残したままの関数を作成すること自体が非常に複雑な文字式の計算となり、非常に時間を要します。たとえば、a=Q+1,b=Q+2、A=a+bならば、A=2*Q+3のような形でAの式を作っているのですが、aとbが非常に複雑で、Aの式の計算に時間を

要するため、先にQ=1とすれば、A=5となることを用いて演算したいのですが、、、

[Kazuya]

なんだか理解が悪くて本当にすいません。まだ何が問題なのかが腑に落ちていません・・。

最初

1. 方程式でQをダイレクトに求める（A(Q) = 0 を解く）は非現実的
2. Qを小数第三位で増分していき、Aがもっとも0に近似した段階のQを求めるという形にしたい。

と伺いました。Q を増分していき、Aを計算する、これは A = f(Q) に該当する計算を、複数の Q について計算することを意味している考えていました。（ちがいますか？）

ただ、その後

• A = f(Q) という Q の値を入力してAを返す関数は、計算コストが高く現実的ではない
• Q=1とすれば、A=5となることを用いて演算したい（？）

とのコメントも頂きまして、混乱しています。

A = f(Q) というのは、最初の質問文内にあるコードからコピペすると

function A = function2getA(Q)
N=2*pi^2;

Z1=(N+2^0.5*Q/2)^0.5;
Z2=(-(-N+2^0.5*Q/2))^0.5;
Z3=(N-2^0.5*Q/2)^0.5;
Z4=(-(-N-2^0.5*Q/2))^0.5;
ah1=(Z1*sinh(Z1)-Z1^2*cosh(Z1))/(2*cosh(Z1)-2-Z1*sinh(Z1));
bh1=(Z1^2-Z1*sinh(Z1))/(2*cosh(Z1)-2-Z1*sinh(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));
ah3=(Z3*sinh(Z3)-Z3^2*cosh(Z3))/(2*cosh(Z3)-2-Z3*sinh(Z3));
bh3=(Z3^2-Z3*sinh(Z3))/(2*cosh(Z3)-2-Z3*sinh(Z3));
a4=(Z4*sin(Z4)-Z4^2*cos(Z4))/(2-2*cos(Z4)-Z4*sin(Z4));
b4=(Z4^2-Z4*sin(Z4))/(2-2*cos(Z4)-Z4*sin(Z4)); 
K1=2*(2+3/2)+(3/4+3*(ah1+a2+ah3+a4)/4);
K2=2*(2+3/2)+(3/4+3*(ah1+a2+ah3+a4)/4);
K3=8+3/2+(ah1+a2+ah3+a4)/2;
K4=-(ah1-a2+ah3-a4)/4;
K5=(ah1-a2-ah3+a4)/(2*2^0.5);
K6=3*(bh1-b2+bh3-b4)/4;
K7=-3/4-(bh1+b2+bh3+b4)/4;
K8=(bh1+b2-bh3-b4)/(2*2^0.5);
K9=(ah1+a2-ah3-a4)/(2*2^0.5);
K10=(bh1-b2-bh3+b4)/(2*2^0.5);
K11=(bh1-b2+bh3-b4)/2; 
K=[K1 K4 K5 K6 K7 K8;K4 K2 K9 K7 K6 K10;K5 K9 K3 K8 K10 K11;...
    K6 K7 K8 K1 K4 K5;K7 K6 K10 K4 K2 K9;K8 K10 K11 K5 K9 K3]; 
A =det(K);

end

となるんだろうと想定していたんですが、違いますか？この計算は文字式処理でもないので、計算自体は重くありません。

[m17td024]

何度も申し訳ございません。

A=det(k)の計算が文字式(Q)を含んだ演算となり、重いと

考えていたのですが、私のパソコンの処理能力が悪いだけなのでしょうか、、、

そこから、solve(A==0,Q)を行ってQの解がでるのがベストなんですが、

いつまで経っても演算が終了しません、、、

[Kazuya]

Q を与えて A の値を計算するのと、文字式処理で A = 0 を解くのはまた全く別の話です。前者は単純な数値計算なので軽いですが、後者は重いでしょうね。。私が今まで提案している手法は前者のみです。

改めて提案としては

1. 単純に Q をいくつか試して A が 0 に近い値になる Q を選ぶ
2. fminbnd関数で試す

どちらも Q を与えて A の値を計算するアプローチなので、計算時間がかかって困るということはないと思いますよ。2 の方がより少ない計算回数で 0 に近くなる Q を算出できます（一般的には）。

[Kazuya]

ここまで来たら・・ということで試してみました。添付の関数をダウンロードして、

Qbest = fminbnd(@function2getAbsA,15,16)

と実行してみてください。一秒もかからず計算は終わると思いますが、結果は

Qbest =
   15.9319

で、その時の A の値（絶対値）は

>> function2getAbsA(15.9319)
ans =
    6.0854

です。0 にはならないみたいなので Q が 15 から 16 までの値をプロットしてみます。

fplot(@function2getAbsA_ver2,[15,16])

となり、この範囲では 0 にはならないみたいですが、一応一番Aが小さくなる Q の値が求まっているみたいですね。（横軸が Q, 縦軸が A の絶対値です）他の範囲でどんな値になるか試してみてください。

function2getAbsA.m と function2getAbsA_ver2 との違いは、、ベクトル入力ができるかどうかです。例えば Q の値を 15から16まで 0.1 刻みで一気に計算できるのは ver2.

例えば function2getAbsA.m でそんな実行をすると、

>> function2getAbsA([15:0.5:16])
エラー:  ^  (line 51)
行列をべき乗にするには次元が正しくありません。行列が正方行列で、べき指数がスカラーであることを確認してくださ
い。行列を要素ごとにべき乗するには、'.^' を使用してください。
エラー: function2getAbsA (line 4)
Z1=(N+2^0.5*Q/2)^0.5; 

とのエラーがでちゃいますが、ver2 なら問題ありません。

>> function2getAbsA_ver2([15:0.5:16])
ans =
   1.0e+06 *
    1.1506    0.4955    0.0708

細かい違いは

https://jp.mathworks.com/help/matlab/matlab_prog/array-vs-matrix-operations.html

を読んでみてください。

[m17td024]

誠にありがとうございます。

非常に精度よく、Aが0に近似する際のQを求めることができました。

お手数お掛けし、申し訳ございませんでした。

[madhan ravi]

+1

[Kazuya]

見つかりましたか！よかったです。