ニュートン法で解が収束しない

こんにちは。
差分式で表した非線形方程式をニュートン法で解いています。が収束しな解あります。ニュートン法は初期値に依存しているため、初期値を可変的にしてみましたがダメでした。何かいい方法はないでしょうか？
参考になるか分かりませんが、使っているプログラムのニュートン法の計算の一部は以下のようです。
call ｇ(x,f,df)
h=f/df
x=x-h
if(dabs(h/x)<1.d-14) then
　return
endif

No.1 回答者： First_Noel 回答日時： 2003/06/11 19:11

これを見ると－１４乗以下で収束ですよね，

変数は倍精度で宣言されていますか？

この回答への補足

さっそくのご回答ありがとうございます。
変数は全て暗黙の倍精度で定義しています。
収束条件を－６乗程度にしても収束しない場合がありした。

補足日時：2003/06/12 01:01

No.2 回答者： noname#108554 回答日時： 2003/06/12 06:12

解が存在するのはまちがいないんですか?

x＾2+1＝0をニュートン法で解くと、いつまでもふらつきますよね。

この回答への補足

ご回答ありがとうございます。
解が存在することはエクセルで確認したので、間違いなくあると思います。

補足日時：2003/06/12 11:42

No.3 回答者： ElectricGamo 回答日時： 2003/06/12 12:07

解いている非線形方程式を書いていただけると一発で答えられます。

まずはNo.2の方がおっしゃってるように解があるかどうかを確認しましょう。そういう時には、収束は遅いですが(最近の計算機は十分速いのでよほど時間のかかる計算を走らせるのではない限り十分だと思いますが)、二分法を使えば便利です。

あるいは解が x=0 だったら収束しないですよ。
その時は、
if ((dabs(x)>=1.d0.AND.dabs(h/x)<1.d-14).OR.
(dabs(x)<1.d0.AND.dabs(h)<1.d-14))
then
　return
endif
とでもしましょう。

この回答への補足

ご回答ありがとうございます。
解いている方程式はとても複雑で、実は連立非線形方程式を差分法を使って解いています。しかし今求めようとしている方程式g(x,f,df)の解が求まれば、他の方程式の解も求まってくる連立方程式になっています。（意味が分からなかったらすみません）そのため連立のニュートン法は使っていません。正確には使いましたがやはり収束しない所がありました。
二分法で確認したところ、解くことができました。
今後の応用も考えてニュートン法でも解けるようにしたいと思っているので、何かいい方法ありましたらアドバイスお願いいたします。

補足日時：2003/06/12 13:01

No.4 回答者： ElectricGamo 回答日時： 2003/06/12 20:23

>二分法で確認したところ、解くことができました。

>今後の応用も考えてニュートン法でも解けるようにし>たいと思っているので、何かいい方法ありましたらア>ドバイスお願いいたします。

具体的な関数の形が分からないと原因に関してやはり正確なことは言えないのですが、私ならば二分法を何回か用いて出したラフな精度の解をニュートン法の初期値として用いて精度を上げます。経験上ほぼ間違いなく収束します。

この回答への補足

再度ご回答ありがとうございます。
二分法で出した解を初期値としてニュートン法を適応してみましたが、やはり収束しませんでした。
もともとこの関数はニュートン法では収束しない性質があるのでしょうか？もしまだ解決策がありましたらよろしくお願いいたします。

補足日時：2003/06/14 15:37

No.5 回答者： ElectricGamo 回答日時： 2003/06/16 04:29

>二分法で出した解を初期値としてニュートン法を適応してみましたが、やはり収束しませんでした。

それはなかなか大変ですねぇ。他にやるとしたら単純にx=x-hとせずにx=x-alpha*hとでもして0<alpha<=1の範囲でfの絶対値が小さくなるようなalphaの一次元探索をしてxを決めていく手があります。

>もともとこの関数はニュートン法では収束しない性質があるのでしょうか？もしまだ解決策がありましたらよろしくお願いいたします。

方程式g(x,f,df)を見てみないことには関数の性質がどうであるとか解決策があるのかどうかは言えません。よろしければここに式を書いてください。

この回答への補足

三度ご回答ありがとうございます。
関数g(x、f、df)は一般的に表すと、
f=x-y*((cos(x)-cos(A))/(x-A))-z*((sin(x)-sin(A))/(x-A))+B
で、少し簡略化しました。
A、Bは定数、ニュートン法で求めたいのはｘです。yとzも本来は変数ですが、ここでは定数扱いです。
何度も何度も申し訳ありません。

補足日時：2003/06/17 11:33

No.6 回答者： ElectricGamo 回答日時： 2003/06/17 19:03

関数の情報ありがとうございます。

試しにこちらでも計算してみますから、収束しない場合の A, B, y, z の具体的な値を一組教えてください。

この回答への補足

お忙しいところ何度もありがとうございます。
関数を少しだけ訂正させていただきます。
f=A*x-y*((cos(x)-cos(B))/(x-B))-z*((sin(x)-sin(B))/(x-B))+C
で、第1項のxに係数をつけただけです。パラメータが、
A=100.0、B=1.22995195D-08、C=-1.24229252D-06
x=1.0、y=7.87169249D-09
の場合収束しませんでした。ですがパラメータが、
A=100.0、B=1.00699949D-08、C=-1.01710309D-06
x=1.0、y=6.44479672D-09
のときは収束しました。
手間をおかけして申し訳ないですが、よろしくお願いいたします。

補足日時：2003/06/18 15:23

No.7 ベストアンサー 回答者： ElectricGamo 回答日時： 2003/06/18 19:04

z の値がないのですが、x → y, y → z と考えてよろしいですよね。

ニュートン法で計算してみると確かに収束しないですね。

与えられたパラメータを入れてグラフを書いてみると解はx = B の近傍にありました。これは図を描かなくても f の両辺を C で割ると分かります。

f/C=A/C*x - y*((cos(x)-cos(B))/(x-B))/C-z*((sin(x)-sin(B))/(x-B))/C + 1

右辺の第二項と第三項はいずれも x = B の近くでゼロに近い値となりますし、A/C * B はほとんど -1です。

で、この関数を見てみると x = B を境にして心電図の波形のように一旦跳ね上がってから潜り込み再び跳ね上がるような形をしています。ちょうど解が存在する辺りに変極点を複数持っています。この場合残念ながらニュートン法ではうまく収束しません。

また、(cos(x)-cos(B))/(x-B)　や (sin(x)-sin(B))/(x-B) は式の上ではx → B の極限で -sin(B) と cos(B) に収束するのですが、数値計算上はゼロで割ることになるので値が不定となります。

解の近くに変極点が複数あること、解の近くで関数の計算自体が困難であることから、この関数に対してニュートン法は適していないと思います。

これだけではあまり前向きのアドバイスにならないので、一言付け加えるとすると、もし解が x = B の近くにあることが前もって分かっているならば、x=Bの近傍で展開した式

f = A*x + y * sin(B) - z * cos(B) + C

を使うという手があります。この場合は手計算で解けちゃいますね。

この回答へのお礼

ご回答ありがとうございます。
ご指摘のとおり、x→y、y→zです。失礼しました。
今回の多くのアドバイスによりニュートン法の性質を理解することができました。
今後の勉強に役立てて行きたいと思います。
長い間本当にありがとうございました。

お礼日時：2003/06/20 12:46