カテゴリー

最新の記事

最近のコメント

最近のトラックバック

月別アーカイブ

ブログ検索

RSSフィード

ブロとも申請フォーム

この人とブロともになる

スポンサーサイト

スポンサー広告
--.--.--
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

ベンチャーブログのランキングに参加しています。
下のバナーをクリックして応援していただけると嬉しいです。
にほんブログ村 ベンチャーブログへ

Newton法で方程式を解く

Scheme/Gauche
2008.12.06
新しい言語のスタディを始めると、やはり簡単なプログラムを作って、これまで使っていた言語との比較を行いたくなる。前回ちょっと言及した東京大学理学部の「Scheme演習」の中に「Newton法で方程式を解け」という問題があったので、やってみることにした。

解を求めたい方程式を

 f(x)=0


とすると、Newton法による方程式の解は以下のように求められる。

Step1.
まず任意のxにおけるf(x)の微分値を求める関数を作成する。この関数を deriv(x) とする。deriv(x)は以下のように定義できる。

 deriv(x) = ( f(x+dx) - f(x) ) / dx



dxが十分に小さい値であれば、deriv(x)の値は十分に微分値に近づくと考えられる。

Step2.
適当なxの値(これを a とする)におけるf(x)の値を求める。つまりf(a)を求めるわけである。

Step3.
x=a における deriv(x) 、つまり deriv(a) を求める。

Step4.
( a - f(a) / deriv(a) ) を求め、これを新しい a とする。

Step5.
Step2~4を繰り返し、f(a)の値を0に近づける。0との差が十分に小さい値になったときに(これを dy とする)、その時の a の値を方程式の解とする。

方程式によっては、最初のaの値によって求められる解が変わってしまったり、解が得られないこともあるが、とりあえず上記の方法で解を求めてみよう。

Scheme/Gaucheでプログラムを作成した場合、コードは次のようになる。

 (define (newton func a dx dy)
  (define (deriv x)
    (/ (- (func (+ x dx))(func x)) dx))
  (cond ((< (func a) dy) a)
    (else (newton func (- a (/ (func a) (deriv a))) dx dy))))


ここで以下のようにf(x)を定義し

 (define (f x)(- (* x x) 2))



以下のようにnewtonを呼び出せば、

 (newton f 5 0.0000000001 0.0000000001)



方程式のひとつの解が得られる。

 gosh> 1.4142135623730887



さてこれを、PHPで書いたらどうなるだろう?
再帰処理を使う場合には、次のようになるはずだ。

 function deriv($f,$x,$dx){
  return((call_user_func($f,$x+$dx)-call_user_func($f,$x))/$dx);
 }

 function newton($f,$a,$dx,$dy ){
  if(call_user_func($f,$a)<$dy){
   return($a);
  }else{
   return(newton($f,$a-(call_user_func($f,$a)/deriv($f,$a,$dx)),$dx,$dy));
  }
 }


さらに f(x) を定義し、newtonを呼び出せばいい。

 function f($x){
  return($x*$x-2);
 }

 echo newton('f',5,0.0000000001,0.0000000001);


Webブラウザからこれを実行すれば、以下のように表示される。

 1.4142135623731



Scheme/Gaucheのように関数内で関数を定義することもできるが、再帰処理を行う場合には内部で定義した関数が重複定義となるため、実行時にエラーになる。そのため微分値を求める deriv 関数は newton 関数の外で定義している。

再帰処理をループに置き換えれば、deriv 関数を newton 関数の中に定義することもできる。

 function newton( $f, $a, $dx, $dy ){
  function deriv( $f, $x, $dx){
   return( ( call_user_func( $f, $x + $dx ) - call_user_func( $f, $x ) ) / $dx );
  }
  while( call_user_func( $f, $a ) >= $dy ){
   $a = $a - ( call_user_func( $f, $a ) / deriv( $f, $a, $dx ) );
  }
  return( $a );
 }


使用法は同じである。

さて、この程度のプログラムであれば、どちらでも簡単に実現できる。ここで私が注目したいのは、どちらの方がシンプルで美しく見えるか、ということである。個人的にはScheme/Gaucheの方がまとまりがいいような気がする。例えばPHPの call_user_func は、ちょっとしたハックの匂いがする。関数名ではなくて、関数そのものを渡せればいいのになあ、と思っちゃうわけだ。

まあこれは、好みの問題なのかもしれませんけどね。
スポンサーサイト

ベンチャーブログのランキングに参加しています。
下のバナーをクリックして応援していただけると嬉しいです。
にほんブログ村 ベンチャーブログへ

FC2Ad

相続 会社設立

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。