ゴールシークのように への返答

投稿で使用できる特殊コードの説明。(別タブで開きます。)
本名は入力しないようにしましょう。
投稿した後で削除するときに使うパスワードです。返答があった後は削除できません。
返答する人が目安にします。相手が小学生か社会人かで返答の仕方も変わります。
最初の投稿が質問の場合、質問者が解決時にチェックしてください。(以降も追加書き込み・返信は可能です。)
※「過去ログ」について書くときはその過去ログのURLも書いてください。

以下の返答は逆順(新しい順)に並んでいます。

投稿者 ひでと  (社会人) 投稿日時 2011/12/20 12:38:28
ぐっさん 様、ありがとうございました。お礼がおそくなりまして、申し訳ありません。

いただいたアドバイスで以下のようにしてみました。変数は省略しています

解を求めたい関数
    Private Function pf(ByVal x As Double) As Double
        Return x ^ 3 + 3 * (pe() - D / 2) * x ^ 2 - 6 * sn * at / b * (pe() + D / 2 - dt) * (D - dt - sn)
    End Function
微分した関数
    Private Function fd(ByVal x As Single) As Double
        Return 3 * x ^ 2 + 6 * (pe() - D / 2) * x
    End Function

解をxnとして
   Private Function xn() As Double
        Dim 初期値 As Double = sd
        Dim dum As Double = newton(10)
        If Math.Abs(pf(dum)) < 0.00001 Then Return dum
        'y=fd*xn+c
        Dim y As Double = pf(dum)
        Dim c As Double = y - fd(dum) * dum
        Dim 初期値1 As Double = (-y - c) / fd(dum)
        Dim 初期値2 As Double = dum
        dum = Bisection(初期値1, 初期値2)
        Return dum
    End Function

まずはニュートン法で
   Private Function newton(ByVal 初期値 As Double) As Double
        'y=fd*xn+c
        Dim xn As Double = 初期値
        Dim i As Integer
        Dim old As Double
        For i = 0 To 100
            Dim y As Double = pf(xn)
            If old = y Then
                Return xn
            End If
            old = y
            If Math.Abs(y) < 0.00001 Then
                Exit For
            End If
            Dim c As Double = y - fd(xn) * xn
            xn = -c / fd(xn)
        Next
        Return xn
    End Function

ニ分法も使ってみて
    Private Function Bisection(ByVal 初期値1 As Double, ByVal 初期値2 As Double) As Double
        If pf(初期値1) * pf(初期値2) > 0 Then
            MsgBox("解を求めることができませんでした")
            Return 初期値2
        End If
        Dim a, b, m As Double
        a = 初期値1 : b = 初期値2
        Dim i As Integer
        For i = 0 To 100
            m = (a + b) / 2
            If Math.Abs(pf(m)) < 0.00001 Then Exit For
            If pf(m) * pf(b) < 0 Then
                a = m
            Else
                b = m
            End If
        Next
        Return m
    End Function

このようにしてみました。
はじめSingleで定義していたとき、ニュートン法だと傾きの精度が足りず、解を求められなかったので、
その近似値?を二分法に使用してみようと考えました。
結局二分法でも精度が足りず、doubleで定義しなおしましたが・・・・・。
試行錯誤、面白かったです。ありがとうございました。
投稿者 ぐっさん  (社会人) 投稿日時 2011/12/2 16:10:24
ニュートン法は、X軸を貫かない極点の付近に初期値を取ると、結果が収束しない場合があります。
少なくとも、f(x)>0となる任意の座標とf(x)<0となる任意の座標を容易に得られるならば、
微分を考えなくてもいい2分法という方法があります。

http://www.math.kobe-u.ac.jp/~taka/asir-book-html/main/node35.html

計算時間は、圧倒的にニュートン法に軍配が上がりますが、人間の感覚からすれば
2分法の計算時間もあっという間です。
投稿者 ひでと  (社会人) 投稿日時 2011/11/29 17:14:04
>shu様 重ねてありがとうございます。
いくつか試行してみます。
数学的なもの自信ないのです・・・。
ただいま49歳ですががんばります。
投稿者 shu  (社会人) 投稿日時 2011/11/29 14:37:15
微分は
簡単に求めるなら(f(x+h) - f(x) )/hを小さいhで求めればよいかと思います。または
すぐれた微分アルゴリズムもあるかと思います。

初期値はこれというのはないので収束性(回数、近い値になっているかなど)を見て適当に変化させる
とよいのではないでしょうか。
投稿者 ひでと  (社会人) 投稿日時 2011/11/29 14:00:57
>shu様ありがとうございます。
ご紹介いただいたURLをちらりと覗いてみました。
この方法だと、関数を微分しないといけないのでしょうか?

微分するのがちょっと辛いので、dを適当な数値にして、xとx+dでの関数のyの値を求めて
その傾きを利用するということでも良いのでしょうか?
そうすると、dの値をどうするか?
初期の値をどうするか?
このあたりのアドバイスなどいただけたらば幸いです。
投稿者 shu  (社会人) 投稿日時 2011/11/29 11:13:17
Googleにてニュートン法
=>
http://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95
http://pc-physics.com/newtonhou1.html
http://akita-nct.jp/yamamoto/lecture/2005/5E/nonlinear_equation/text/html/node4.html

など
投稿者 ひでと  (社会人) 投稿日時 2011/11/29 09:21:49
変数xに関する関数f(x)があります(式は分っています)。
f(x)=0となる解を求める簡単な方法を教えてください。
Excelのゴールシークのようなイメージです。
建築の柱脚部の強度計算に必要なものです。
できれば VB2005での具体的なものがありがたいです。