高等学校数学B/数値計算とコンピューター

出典: フリー教科書『ウィキブックス(Wikibooks)』
ナビゲーションに移動 検索に移動
高等学校数学 > 高等学校数学B > 高等学校数学B/数値計算とコンピューター

数値計算とコンピューター[編集]

初等的な算法を扱い、計算機を用いてそれを計算する方法を学ぶ。プログラム例としては、PythonSchemeという言語を扱うが言語の詳細に立ち入らず、考え方を学ぶことが重要となる。

整数の算法[編集]

ユークリッドの互除法[編集]

ユークリッドの互除法は2つの整数の最大公約数を求める算法である。ある整数m, n (m > n > 0) をとる。このときユークリッドの互除法は

  1. mをnで割った余りを計算し、それをrとおく。このときr=0なら3に進み、r 0なら、2に進む。
  2. mを以前のnの値で置き換え、nをrの値で置き換え、1に戻る。
  3. nの値が最大公約数となっている。

で与えられる。

(導出)

m,nが互いに素であるときを考える。mをnで割った商をa、余りをrとするとき、

(rnm)

が成り立つ。ここで、仮にn、rが共通因数を持つならその因数はmの因数でもあるがこれはm、nが互いに素であることに矛盾する。よって、n、rは互いに素である。ここから上の1、2を行なうと互いに素でありより小さい2つの整数n,rが得られる。これを繰りかえすと小さい側の整数は1となる。

実際余りが2以上になるときは2数が互いに素であることから、次の計算で更に小さい数が得られ、余りが0になることは小さい方の数が1である場合を除いて、2数が互いに素であることに反する。よって、確かに小さい側の整数は1となる。よって、m,nが互いに素であるときユークリッドの互除法は確かめられた。次にm,nが最大公約数Mを持つときを考える。このときもmをnで割った商をa、余りをrとするとき、

(rnm)

が成り立つが、

を考えると、rもm、nと同じ最大公約数Mを持つ。r,m,nをMで割ったものをそれぞれr',m',n'とおくと、これらは互いに素であるが(最大公約数の定義)、このとき上の2数が互いに素であるときのユークリッド互除法の導出から小さい方の整数は1が得られる。よって元の整数に戻るためにMをかけることで、この方法が2数の最大公約数Mを与えることが分かる。よって、m、nが共通因数を持つ場合にもユークリッド互除法は示された。

実際の計算には計算機を用いると(特に2数が大きいときには)便利である。

Pythonによるプログラム例

def euclid():
    left =45 
    right = 30
    assert (left > 0 and right > 0),"not a positive number"
    def exactly_divides_the_other():
	if (left > right):
	    if (left % right == 0):
		return True
	elif (left < right):
	    if (right % left == 0):
		return True
	return False 
    while(not exactly_divides_the_other()): 
	if (left > right):
	    left = left % right
	elif (right > left):
	    right = right % left
    #Then we get a mcd. The smaller number is that.
    if (left < right):
	return left
    elif (left > right):
	return right

print euclid()
#ok that work well if I put 
#(left ,right)= (45 ,30)
#(45 ,28),(30 ,28)

Schemeによるプログラム例

(define (euclid m n)
  (let ((r (modulo m n)))    
    (if (zero? r)             ;ここまでが導出過程の1
        n                     ;ここが導出過程の3
        (euclid n r))))      ;ここが導出過程の2

;;;実行例
;;> (euclid 45 30)
;;15
;;> (euclid 45 28)
;;1
;;> (euclid 30 28)
;;2

実数の算法[編集]

2分法[編集]

ある関数f(x)とx軸との接点を求める方法の1つに、2分法がある。特にf(x)が求める点で正の傾きを持っているものとして考える。この方法は、

  1. 範囲[a,b]内にx軸と求める関数f(x)の接点が含まれるように、2数a,bを定める。
  2. mid_point = (a+b)/2 を計算する。もしf(mid_point)が十分に0に近ければ4に進む。
  3. もしf(mid_point)0なら、mid_pointの値をbの値に代入し、2に戻る。もし、f(mid_point)0なら、mid_pointの値をaの値に代入し、2に戻る。
  4. mid_pointの値が求める接点のx座標である。

この方法は元々の範囲[a,b]の中点を取り、解が中点から見てどちらにあるかを判断し、範囲を狭めていく方法である。

Pythonによるプログラム例

def bisection():
#Could be a too simple function. But seems ok,
#it just is an example.
    def func(x):
	return x**2-1
#2 magic variables
    left = 0.0
    right  = 3.0
    assert (type(left) == float
	    and type(right) == float),"That is not real numbers"
    mid = (left +right)/2
    while not (-1.0e-10 <func(mid) < -1.0e-10):
	mid = (right + left)/2
	if (func(mid) < 0):
	    left = mid
	elif (func(mid) > 0 ):
	    right = mid
#That could hardly be true ... .
	else:
	    return mid
    return mid

print bisection()


#ok that work well if I put 
#func (x):
#    x-1 or x**2 - 1
#The result of both calculation was 1.0.


このコードはf(x)=x-1、または、f(x)=のときに試された。結果はどちらもx=1.0となり、正しい値を返している。

Schemeによるコード例

(define (bisection f a b)                       ;手順1。
  (let ((e (expt 10 -10))               
        (mid_point (/ (+ a b) 2)))             ;手順2。中点の計算。
    (cond ((or (zero? (f mid_point)) 
               (< (- e) (f mid_point) e)) 
           (exact->inexact mid_point))         ;ここまでが手順4。
          ((> (f mid_point) 0) 
           (bisection f a mid_point))
          (else (bisection f mid_point b)))))  ;ここまでが手順3

;;;実行例
;;> (bisection (lambda (x) (- x 1)) 0 3)           ;x-1の解を0〜3間で探す。
;;0.9999999999417923
;;> (bisection (lambda (x) (- (expt x 2) 1)) 0 3)  ;x^2-1の解を0〜3間で探す。
;;1.0000000000291038

このコードもλ(x)=x-1、または、λ(x)=のときに試された。Python版と違うが、結果はどちらもx=1.0に極めて近い値を返している。計算精度の違い、である。

台形公式[編集]

台形公式は、あるグラフf(x)とx軸とx=a,x=bに囲まれた面積を近似的に求める公式である。この公式では、[a,b]の範囲をN個の小さい範囲に分け、i個目の範囲を、と書く。このときその範囲においては求める面積を台形で近似しても面積のずれは小さい。

正確な面積と台形の面積のずれの絵

ここで、台形の面積

で書かれることを考慮すると、求める面積Sは、

で近似できることが分かる。

Pythonによるプログラム例では、半径1の円の面積を近似的に求め、それによっての値を計算する。

from math import sqrt,pi

def trapezoid_formula():
    def func(x):
	return sqrt(1-x**2)
    sum = 0.0
    a = 0.0
    b = 1.0
    assert (type (sum) ==float and
	    type (a) ==float and
	    type (b) ==float), "not a real number" 
    N = 20
    dx = (b-a)/N
    for i in range(N):
#Section of trapezoid ... .
	sum += (func(a+dx * i)+ func(a+dx*(i+1)) ) * dx /2
    return sum

print trapezoid_formula()
print pi/4

#ok that work well
#the result reads 
#0.782116219939 for  trapezoid_formula()
#0.785398163397 for pi/4 


実際のの値と近い値が得られていることが分かる。

Schemeによるプログラム例

(define (trapezoid_formula f a b)
  (let ((n 20))
    (let ((dx (/ (- b a) n)))
      (let loop ((i 0) (sum 0))
        (if (= i n)
            (exact->inexact sum)
            (loop (+ i 1) 
                  (+ sum (* (+ (f (+ a (* dx i)))
                               (f (+ a (* dx (+ i 1)))))
                            (/ dx 2)))))))))

;;;実行例
;;> (trapezoid_formula (lambda (x)
;;                       (sqrt (- 1 (expt x 2))))
;;                     0 1)
;;0.7821162199387454
;;> (/ pi 4)
;;0.7853981633974483

こちらも実際のの値と近い値が得られていることが分かる。