×[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
prev>>Lispの真偽判定｜ HOME ｜next>>2015/01/11/ミドリムシで、空を飛ぶ。
指数関数
xを任意の非負実数、mを任意の非負整数とするとき、xのm乗の再帰関数mypwrを定義します。
if関数
****
(defun mypwr (x m)
(if (= m 0)
1.0
(* x (mypwr x (- m 1)))))
; [1]> (mypwr 2 3)
; 8.0
; [2]> (mypwr 0.5 3)
; 0.125
****
cond関数
****
(defun mypwr (x m)
(cond ((= m 0) 1.0)
(T (* x (mypwr x (- m 1))))))
; [1]> (mypwr 2 3)
; 8.0
; [2]> (mypwr 0.5 3)
; 0.125
****
xが2、mが3の場合について考えてみましょう。if文で"(= m 0)"と記述されているので、条件に当てはまらないときの手続きである"(* x (mypwr x (- m 1)))))"を行います。"(* 2 (mypwr 2 2))))"と解釈ができるので"(mypwr 2 2)"について考えます。同様にmは0ではないので"(* 2 (mypwr 2 1))))"となり、"(mypwr 2 1)"でもmは0ではないので"(* 2 (mypwr 2 0))))"となります。"(mypwr 2 0)"ではmの値が0なので1.0を返し、"(mypwr 2 1)"は1.0なので置き換えると"(* 2 1.0)"となります。また同様に繰り返していくと"(mypwr 2 2.0)"、"(mypwr 2 4.0)"となり最終的な答えは8.0になります。
平方根
任意の正実数cの平方根c^(1/2)を求める関数mysqrtを定義してみましょう。実数xを引数とする実関数をF(x)=x^(2)-cとするとき、F(x)が0ならば等式の関係が成り立ちます。xとcをそれぞれ2乗している理由は最終的な答えをxでだすためです。x^(2)=cならばxを返せばその値はc^(1/2)になります。さて、この方程式を求めるにはニュートン法と呼ばれる再帰的な計算方法でその近似値解を求めます。初期値を正の実数x(0)にしたときの近似値解x(n)は
x(n)=x(n-1)-F(x(n-1))/f(x(n-1))
になります。x(0)はxの最初の解、x(n)は上の式を再帰的にn回繰り返したときに求められる解のことです。その2つの解の間にはx(1),x(2),...,x(n-2),x(n-1)とあります。さらに、この式でのf(x(n-1))はF(x(n-1))の導関数を表します。つまり、f(x(n-1))はF(x(n-1))を微分したものになります。これらのことをF(x)=x^(2)-cに置き換えてみると
x(n)=x(n-1)-(x(n-1)^(2)-c)/(2x(n-1))
になり、分数部分を分解して再度整理すると
x(n)=x(n-1)/2+c/(2x(n-1))
Lispで記述すると
****
(defun better (x c)
(+ (/ x 2.0)
(/ c (* 2.0 x))))
****
になりました。さて、ニュートン法とはnの値が大きくなればなるほど解に収束していきますが計算を永遠に続けていくわけにはいきません。どこかで誤差(妥協)を決めなければならないのです。それは人それぞれですが、今回は誤差0.001で考えてみましょう。
実際の解F(x)と再帰的に求めた解F(x(n))の誤差を求めるにはその２つの解の絶対値を求めて、それが誤差未満かどうかを判断します。絶対値を求める関数は前回の真偽の判定で定義した関数absoluteを使用します。
****
(defun absolute (x)
(cond ((< x 0) (- 0 x))
(T x)))
****
関数absoluteで求めた値が誤差より小さければ計算を止めるように命令をすればよいので、ここでは"T"を返すようにします。誤差を変数errと置くと
|F(x)-F(x(n))|<err
F(x)を変数a、F(x(n))を変数bと置いてLispで記述したときの関数nearは
****
(defun near (a b err)
(< (absolute (- a b)) err))
****
になります。最後に関数nearで"T"が帰ってきた場合はxを返し、"NIL"が帰ってきた場合はxを再度計算します。つまり、xの値を関数betterでさらに計算した値に置き換える処理をします。上記の処理をする関数の名前をmysqrt-coreとすると
****
(defun mysqrt-core (c x err)
(if (near c (* x x) err)
x
(mysqrt-core c (better x c) err)))
****
になり、関数mysqrtは
****
(defun mysqrt (c)
(mysqrt-core c 1.0 0.001))
****
と記述ができます。
****
(defun better (x c)
(+ (/ x 2.0)
(/ c (* 2.0 x))))
(defun absolute (x)
(cond ((< x 0) (- 0 x))
(T x)))
(defun near (a b err)
(< (absolute (- a b)) err))
(defun mysqrt-core (c x err)
(if (near c (* x x) err)
x
(mysqrt-core c (better x c) err)))
(defun mysqrt (c)
(mysqrt-core c 1.0 0.001))
; [1]> (mysqrt 2)
; 1.4142157
; [2]> (mysqrt 5)
; 2.2360687
; [3]> (mysqrt 4)
; 2.0
演習
任意の正実数cと2以上の任意の整数mが仮引数として与えられたときのcのm乗根(>0)を計算する関数rootを定義し適当な値で評価せよ。なお、初期値xは1.0、誤差errは0.001とし、必要な関数は自分で定義をすること。
演習回答
****
(defun better (x c m)
(+ (* (- 1 (/ 1 m)) x)
(/ c (* m (mypwr x (- m 1))))))
(defun absolute (x)
(cond ((< x 0) (- 0 x))
(T x)))
(defun near (a b err)
(< (absolute (- a b)) err))
(defun mypwr (x m)
(if (= m 0)
1.0
(* x (mypwr x (- m 1)))))
(defun root-core (c m x err)
(if (near c (mypwr x m) err)
x
(root-core c m (better x c m) err)))
(defun root (c m)
(root-core c m 1.0 0.001))
; [1]> (root 3 8)
; 1.1472127
; [2]> (root 4 2)
; 2.0
; [3]> (root 15 4)
; 1.9679902
演習解説
この問題でポイントになるのはF(x)からx(n)が求められるか、x(n)をLispで記述できるかの2点です。それさえできれば2乗根と同じ要領でできます。
まず、cのm乗根をF(x)=x^(m)-cとするときf(x)はm*x^(m-1)となります。そしてx(n)はニュートン法の式に代入していきx(n)=x(n-1)-(x(n-1)^(m)-c)/(m*x(n-1)^(m-1))を得ます。式を整理してみるとx(n)=(1-(1/m))*x(n-1)+c/(m*x(n-1)^(m-1))となりx(n)を求めることができました。
次にx(n)をLispで記述するには、xのm乗を計算しなくてはなりません。指数の計算に少し戸惑った人がいると思いますがこのページの最初の定義式に指数を計算する関数mypwrがあるのでそれを利用します。
[0回]
PR
prev>>Lispの真偽判定｜ HOME ｜next>>2015/01/11/ミドリムシで、空を飛ぶ。
