反復計算で指数方程式の解を求めたい

よろしくお願いいたします。以下方程式（単調増加関数）のｘについて求めたいです。

a x + (b x)^1/(n+1) +(c x)^1/(m+1) - d=0

a, b, c, d （定数）を各列に代入、各行の右側にxを求める式を代入し計算するシートを作りたいです。反復計算、循環参照を使用しa,b,c,dの列の横に数式をコピペして一括て求めたいです。
近似解でも良いです。

ためしに式を変形し
x={(b x)^1/(n+1) +(c x)^1/(m+1) - d}/a
として、たとえば
a,b,c,dがA1, B1, C1, D1に代入されており

E１：=((B1*F1)^1/(n+1)+(C1*F1)^1/(m+1) -D1)/A1
F1：=E1

として循環参照をさせ、反復回数100、精度0.0001に設定し計算したところ
どちらのセルも計算されず、＃NUM となってしまいました。
（ｄ＝０のときのみｘ＝０と計算できました。）

a～d及び指数のn,mの数値の範囲は以下です。
a：0.04-0.07
b：90-100
c：70-90
d：0-100
n：:0.157
m：0.352

ニュートン法を用いると負のべき乗または負の指数を計算することになって＃NUMとなってうまくいきません。この問題をどうにか回避して計算する方法をご存知の方、ご教示いただきたくお願いします。

No.6 ベストアンサー 回答者： nag0720 回答日時： 2013/05/01 01:05

＞ａ，ｂ，ｃ，ｄの各条件が1000行程あり、各条件毎にｘを算出するには、

No.5で書いたH1,I1,G2を縦ではなく横に並べればいい。

３つ並べるのが面倒なら、１つの式にまとめてください。

A～Fにa,b,c,d,n,m
Gにxの初期値
H1に、=G1-($A1*G1+($B1*G1)^(1/($E1+1))+($C1*G1)^(1/($F1+1))-$D1)/($A1+$B1/($E1+1)*($B1*G1)^(1/($E1+1)-1)+$C1/($F1+1)*($C1*G1)^(1/($F1+1)-1))

I,J,K,L,M,N,・・・・にHをコピー

10回くらいもやれば収束するでしょう。

この回答へのお礼

早速のご回答ありがとうございます。
ばっちり計算できました。
大変助かりました。
ありがとうございました。

お礼日時：2013/05/01 02:08

No.5 回答者： nag0720 回答日時： 2013/04/30 22:18

ニュートン法です。

(a) A1：0.06294
(b) B1：92.674
(c) C1：80.531
(d) D1：0.45
(n) E1：0.157
(m) F1：0.352
(x0) G1：0.0001　（初期値）
(f) H1：=A1*G1+(B1*G1)^(1/(E1+1))+(C1*G1)^(1/(F1+1))-D1
(f') I1：=A1+B1/(E1+1)*(B1*G1)^(1/(E1+1)-1)+C1/(F1+1)*(C1*G1)^(1/(F1+1)-1)

(x1) G2：=G1-H1/I1

他のセルは、上のセルをコピー。

これで６回目ぐらいで、x=0.001787631　に収束します。


この関数は上に凸の増加関数なので、初期値の取り方に注意する必要があります。
大きい値を初期値にするとxがマイナスになってしまいます。
できるだけ0に近い正値を初期値にしたほうがいいでしょう。

この回答へのお礼

ご回答頂きありがとうございます。
ご指摘の通り、初期値の取り方次第で変わってくることを理解できました。
ありがとうございます。

ａ，ｂ，ｃ，ｄの各条件が1000行程あり、各条件毎にｘを算出するには、どうしたらよいかご存じありませんでしょうか。
（1000回同じことをすればよいのでしょうが、、、、。）
ご教示いただければ幸いです。

お礼日時：2013/04/30 22:56

No.4 回答者： spring135 回答日時： 2013/04/30 21:48

＃3です。

プログラムにミスがあって結果がおかしくなっていました。微分係数に不安定があるようで収束が悪いですが、とりあえず1000回繰り返したところ収束値だけは出ました。

x=0.285632676376962

この回答へのお礼

ありがとうございます。
できれば、各ａ，ｂ，ｃ，ｄのシートに対しxを一括計算させるようにしたいです。
計算数が多く（1000行くらいあり、また多くのシートを処理しなければなりません。）
ａ，ｂ，ｃ，ｄが各列に入力されているシートに各行ごとにｘを計算するにはどうしたらよいか、
ご存じであればご教示ください。

お礼日時：2013/04/30 22:41

No.3 回答者： spring135 回答日時： 2013/04/30 21:39

まずはy=a x + (b x)^1/(n+1) +(c x)^1/(m+1) - dのグラフを書いてみてください。

そしたら一目瞭然です。


yを微分してy'を求め

ニュートン－ラフソンで

x(n)=x(n-1)-y/y'

で簡単に求めることができます。


y=a x + (b x)^1/(n+1) +(c x)^1/(m+1) - dがほとんど直線状なので一発で収束します。


多分(b x)^1/(n+1) +(c x)^1/(m+1)は効いてこないのでしょう。

a：0.05
b：95
c：80
d：50
n：:0.157
m：0.352

のとき

x=-0.315623356


＞ためしに式を変形し
x={(b x)^1/(n+1) +(c x)^1/(m+1) - d}/a


何のための変形ですか。

この回答へのお礼

ご回答をありがとうございます。
おっしゃる通り確かにグラフを描くと一目瞭然なのですが、
なぜか収束できません。

x={(b x)^1/(n+1) +(c x)^1/(m+1) - d}/a

の変形は、循環参照させてセル同士を参照させると同じ値に収束して求めることができた経験からです。

エクセルの循環参照と反復計算のヘルプを見て可能と思い実施してみたところうまく収束しました。
その時は、指数が整数のためか収束し解を得ることができました。

今回は、ニュートン法の途中でべき乗の中身が負となる場合がありエクセルが計算できず
＃NUMを返したとようです。
（負のべきはエクセルでサポートされていないようです。）

お礼日時：2013/04/30 22:39

No.2 回答者： ddtddtddt 回答日時： 2013/04/30 18:42

＞ニュートン法を用いると・・・

と書かれているので、ゴールシークだけでなく、マクロ（VBA）もやった事があると思って良いですか？。以下は、VBAを使えるが前提です。

　単調増加になるa，b，c，d，n，mを与える訳ですから、解はあったとして1個です。2分法（bysection）はどうでしょう？。解が0以上にあるなら、計算過程で負数のべき乗や負数の指数を計算する事はありません。

　2分法はニュートン方に比べて遅いと言われていますが、最近のPCは信じられないくらい速いので、この程度なら恐らく、体感速度は0ですよ。それにアルゴリズムは簡単だし（与式を単に計算するだけ(^^)）、解精度の制御もしやすいし・・・(^^)。

この回答への補足

申し訳ありません。
2分法を良く知らず、勝手なコメントをしてしまいました。
（引数に関数と、ｘを入れてｘの値を返すVBAなるものをお尋ねしてしまいました。）
2分法の計算方法とVBAをネットで探してみましたが、なかなか奥深く見つけるのが大変でした。
勝手なコメントをして申し訳ありません。
もし、お気を悪くされず、質問の対策方法を引き続きご教示いただければ幸いです。

補足日時：2013/04/30 19:59

この回答へのお礼

ご回答頂きありがとうございます。
おっしゃるとおり、単調増加ですので試算すると解は1個なのですが、質問の結果となってしまいました。
2分法については、知見がないのです。ご教示いただければ幸いです。
当方PC環境はWIN7でOFFICE２０１０です。
VBA等で関数を定義できて、引数に方程式とｘを入れて、解くマクロとかあると助かります。
ご存じであれば、自分はVBAの知識はほとんどないのですが、ご教示いただけると助かります。

お礼日時：2013/04/30 19:47

No.1 回答者： nag0720 回答日時： 2013/04/30 16:45

a x + (b x)^(1/(n+1)) +(c x)^(1/(m+1)) - d = 0　でしょうか。

＞ニュートン法を用いると負のべき乗または負の指数を計算することになって

xが負になるのならだめですが、指数は負でも計算できます。

ニュートン法で問題なく計算できますよ。

a,b,c,dがどんな値のとき＃NUMになったんでしょうか。

この回答へのお礼

早速のご回答ありがとうございます。
ニュートン法でためしにx_k=x_k-1-f(x)/f'(x)で計算しましたが、初期値を１にしたらいきなり、
xが負で帰ってきてしまいました。以降f'(x)が＃NUMとなって以降の計算ができませんでした。
値は質問の範囲で求めたいですが、ためしに以下条件で行いました。

a：0.06294
b：92.674
c：80.531
d：0.45
n：:0.157
m：0.352

お礼日時：2013/04/30 19:41