私的数学日誌

17日目「Newton法についてとそのプログラム」

皆さんはNewton法という数値解を見つける方法を知っているだろうか。

非常に簡単な考え方をつかって数値解を求めるが、解への収束が結構早いことが特徴の数値解析の一種の方法である。とりあえず、考え方を説明しよう。

とりあえず、解を求めたい関数を図に書いてみる。今回のプログラムでは、√を求めたいので、f(x)=x^2-mとでもおいておこう。ここで、ある点を初期値としてとる。この初期値は、ある程度解に近い必要がある。なぜならば、もう一つ解があったとしたら、そちらの方に収束する可能性もあるからである。また、そもそも収束をしない点もあるので注意されたい。しかし、f(x)では収束しないということはほとんど起こりえないので、初期値として適当にmをとっておく。x_0=mとおいてしまおう。

これで準備は完了した。図に書いた関数に初期値をプロットする。そして、そこから接線を引く。すると、その接線のx軸との交点は初期値よりも求めたい解に近づいていることが、図より明らかにわかる。

このことを利用して、接線の方程式を用いて、数列の漸化式をつくり、極限が解に収束していくようにしていく。

接線の方程式はx=aを接点とするとき、g(x)=f'(a)(x-a)+f(a)とかかれ、x_nを接点としてx_(n+1)を求めていくので、f'(x_n)(x_(n+1)-x_n)+f(x_n)=0より、x_(n+1)={f'(x_n)x_n-f(x_n)}/f'(x_n)と漸化式が立てられる。x_(n+1)=x_n-f(x_n)/f'(x_n)

とできて、f'(x)=2xより、x_(n+1)=1/2(x_n+m/x_n)とできる。x_0=m。

この数列はかなり早く解に近づいていく。pythonでテキトーに組んでみたプログラムがこれだ。

u=input('input number')

t=int(u)

def nf(k):

if k==0:

return t

else:

return 0.5*(nf(k-1)+t/nf(k-1))

for i in range(0,15):

print(i,nf(i))

例えば、2をこれに入れてみる。すると5番目からは限界の制度に到達していることがわかる。数字が大きいほど、十分な精度が得られるまでは時間がかかるが、その上昇量も少なめである。

ちゃんとプログラムを組むのなら、終了の条件を、隣り合った項の差によって決める方がいいと思う。あと、python初心者の私でもできたので、練習用はよいと思う。