旧課程(-2012年度)高等学校数学B/数値計算とコンピューター

高等学校数学 > 高等学校数学B > 旧課程(-2012年度)高等学校数学B/数値計算とコンピューター

数値計算とコンピューター

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

整数の算法

ユークリッドの互除法

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

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

で与えられる。

(導出)

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

$m=na+r$ ただし $(r<n<m)$

が成り立つ。ここで、仮に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とするとき、

$m=na+r$ ただし $(r<n<m)$

が成り立つが、

$r=m-na$

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

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

Pythonによるプログラム例

def euclid(m, n):
    print(f"euclid({m}, {n})")
    if (n == 0):
        return m
    return euclid(n, m % n)

print(euclid(45,30))
print(euclid(45,28))
print(euclid(30,28))

実行結果

euclid(45, 30)
euclid(30, 15)
euclid(15, 0)
15
euclid(45, 28)
euclid(28, 17)
euclid(17, 11)
euclid(11, 6)
euclid(6, 5)
euclid(5, 1)
euclid(1, 0)
1
euclid(30, 28)
euclid(28, 2)
euclid(2, 0) 
2

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によるコード例

from math import isfinite

def bisection(func, left: float, right: float) -> float:
    # acceptance inspection
    assert (callable(func)),"func is not callable."
    assert (isfinite(left)),f"The left({left}) is not a finite number."
    assert (isfinite(right)),f"The right({right}) side is not a finite number."
    assert (left <= right),f"The left({left}) is bigger than the right({right})."
    
    # Implementation of core algorithms
    def core(f, low: float, high: float) -> float:
        x = (low + high) / 2
        fx = f(x)
        if (abs(fx) < +1.0e-10):
            return x
        if fx < 0.0:
            low = x
        else:
            high = x
        return core(f, low, high)
        
    return core(func, left, right)

print(bisection(lambda x: x-1, 0.0, 3.0))
print(bisection(lambda x: x*x-1, 0, 3))

実行結果

0.9999999999417923 
1.0000000000291038

このコードは$\lambda (x)=x-1$、または、$\lambda (x)=x^{2}-1$のときに試された。結果は 0.9999999999417923 および 1.0000000000291038 であり、充分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

(print (bisection (lambda (x) (- x 1)) 0 3))           ;x-1の解を0〜3間で探す。
(print (bisection (lambda (x) (- (expt x 2) 1)) 0 3))  ;x^2-1の解を0〜3間で探す。

実行結果

0.9999999999417923 
1.0000000000291038

このコードも$\lambda (x)=x-1$、または、$\lambda (x)=x^{2}-1$のときに試された。Python版の結果と一致している。

台形公式

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

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

ここで、台形の面積$s_{i}$は

$s_{i}={\frac {1}{2}}\{f(x_{i})+f(x_{i+1})\}\cdot (x_{i+1}-x_{i})$

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

$S=\sum _{i=0}^{N}s_{i}$

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

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

trapezoid.py

from math import sqrt,pi
from numbers import Number

def trapezoid_formula(f, a, b):
    assert callable(f), "f must be a callable"
    assert isinstance(a, Number), "a must be a number"
    assert isinstance(b, Number), "b must be a number"
    
    n = 20
    dx = (b - a) / n
    sum = 0
    for i in range(n):
        sum += (f(a + dx * i) + f(a + dx * (i + 1))) * dx / 2
    return sum

print(trapezoid_formula(lambda x: sqrt(1 - x ** 2), 0, 1))
print(pi/4)

実行結果

0.7821162199387454
0.7853981633974483

実際の$\pi /4$の値と近い値が得られていることが分かる。

Schemeによるプログラム例

trapezoid.scm

(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)))))))))
                        
(print (trapezoid_formula (lambda (x)
                       (sqrt (- 1 (expt x 2))))
                     0 1) )

(print (atan 1.0))

実行結果

0.7821162199387455 
0.7853981633974483

こちらも実際の$\pi /4$の値と近い値が得られていることが分かる。