Matlabについて質問2

二分法とニュートン法について、どのようにプログラムを書けばいいのでしょうか？
具体的には、fortranでは、以下のように書かれているのですが、、、
ｙ＝-x^2+5x+2

（以下のページ
http://www8.plala.or.jp/uemura/sim2.html
を参照しています。）

（二分法）
DOUBLE PRECISION X,AX,X1,X2,Y1,Y2,YA,Y
X=1
H=0.1
X1=0　!低い方のxの初期値
X2=5　!高い方のxの初期値
100 AX = (X1 + X2) / 2　!２つの初期値の平均
Y1 = Y(X1)
Y2 = Y(X2)
YA = Y(AX)
IF (YA.EQ.Y1.AND.YA.EQ.Y2) GOTO 200
IF (Y1.GE.Y2) THEN
　X2 = AX
　GOTO 100
END IF
IF (Y2.GE.Y1) THEN
　X1 = AX
　GOTO 100
END IF
200 WRITE(*,*) 'X*= ', X1
WRITE(*,*) 'Y*= ', Y1
END

DOUBLE PRECISION FUNCTION Y(X)
DOUBLE PRECISION X
Y=-X**2+5*X+2
END

これを実行すると以下のようになる。

X*= 2.499999981373549
Y*= 8.250000000000000

（ニュートン法）
DOUBLE PRECISION X,H,DY1,DY2,DY3,Y1,Y,DY
X=1　!xの初期値
H=0.1　!ステップ幅の初期値
100 DY1 = DY(X)
DY2 = DY(X + H)
DY3 = DY(X - H)
IF (ABS(DY1).LT..000001) GOTO 200　!収束判定
IF (ABS(DY2).LT.ABS(DY1).AND.ABS(DY2).LT.ABS(DY3)) THEN
　X = X + H
　GOTO 100
END IF
IF (ABS(DY3).LT.ABS(DY1).AND.ABS(DY3).LT.ABS(DY2)) THEN
　X = X - H
　GOTO 100
END IF
IF (ABS(DY1).LE.ABS(DY2).AND.ABS(DY1).LE.ABS(DY3)) THEN
　H = H / 10
　GOTO 100
END IF
200 Y1 = Y(X)
WRITE(*,*) 'X*= ',X
WRITE(*,*) 'Y*= ',Y1
END

DOUBLE PRECISION FUNCTION Y(X)
DOUBLE PRECISION X
Y=-X**2+5*X+2
END

DOUBLE PRECISION FUNCTION DY(X)　!微分式の定義
DOUBLE PRECISION X
DY=-2*X+5
END

これを実行すれば以下のようになる。

X*= 2.500000022351742
Y*= 8.250000000000000

No.1 回答者： Kules 回答日時： 2012/01/06 00:18

なんだかなあ…

これってMatlabの問題なんですかね？
二分法だけ取り上げます。
あなたの参照したページの中では、
「解x*よりも小さい値と大きい値の２つの初期値を与え、その平均値を計算し、この３つのxから得られるyの中で最も小さいxを捨て、新たに平均値を計算するという技法である」
って書いてますね？
また、かのwikipediaでは
http://ja.wikipedia.org/wiki/%E4%BA%8C%E5%88%86% …
ここでは、f(x) = 0となるxを求める方法について説明する。
1.f(x1)とf(x2)とで符号が異なるような区間下限x1と区間上限x2を定める。
2.x1とx2の中間点xMを求める。
3.f(xM)の符号がf(x1)と同じであればx1をxMで置き換え、f(x2)と同じであればx2をxMで置き換える。
4.2.に戻って操作を繰り返すことにより、f(x) = 0となるxに近づく。
という説明が書かれています。

（Matlab独特の書き方に頼るということをするのでなければ）この文章を愚直にプログラムにおこすだけの話だと思うんですが。プログラムとして難しそうなのは上記手順の3.でしょうか。これもif文の基本的な構文を知っていればできますね。
ちなみに、Matlabには入力値の符号を返す関数signがあるので、これを有効活用してやればよいでしょう。

ニュートン法も同じです。「文章で書かれたことをマジメにプログラムにおこす」この一手です。

「とりあえずこう書いてみたんだけど動かない」とか「こう書いたんだけど時間がかかってしょうがない、もっと高速に出来ないものか」とかであればその点補足いただければより有益な回答が付くのではないかと思われます。

途中厳しいことも書きましたが、参考になれば幸いです。