LK法のパフォーマンス

こちらの記事では、1次元の場合についてLK法を実装し、パフォーマンスについて解説します。

この技術についての投稿は、一つの索引用の記事にまとめられています。全体を俯瞰する場合はそちらを御覧ください。

m-yoshi-1700.hatenablog.com

前々回からLK法の式の導出をしてきました。

ここまでで求めた式は以下になります。

${ \displaystyle \label{eq10} h_{0} = 0\\ h_{k + 1} = h_{k} + \left[ \frac {\sum_{x} w(x) F'(x + h_{k})[G(x) - F(x + h_{k})]} {\sum_{x} w(x) F'(x + h_{k})^{2}} \right] \tag{10} }$

こちらの式で実際のhに収束する条件と、収束するための速度について考えます。

LK法の実装

ここで、F(x)とG(x)を以下の様に置きます。

${ \displaystyle F(x) = \sin x\\ G(x) = F(x + h) = sin(x + h) }$

こちらをpythonで実装します。まず、F(x)とG(x)をグラフに描画してみます。

#!/usr/bin/env python
#-*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

# hの真値
h_g = 0.3


def f(x):
    return np.sin(x)


def g(x):

    return np.sin(x + h_g)


def main():
    x = np.arange(0, 2 * np.pi, 0.01)

    f_val = f(x)

    g_val = g(x)

    plt.plot(x, f_val, label="f(x)")
    plt.plot(x, g_val, label="g(x)")
    plt.legend()
    plt.show()

if __name__ == "__main__":
    main()

こちらを実行すると、次のようなグラフが表示されます。

f:id:m-yoshi-1700:20170319122604p:plain

Fが正弦波で、Gがそれをhだけずらしたものになっています。

ここから、式\eqref{eq10}を用いてhを求めます。

def weight(x):
    diff = np.absolute(g_1st_deriv(x) - f_1st_deriv(x))
    if diff == 0:
        return 0
    return (1 / diff)

ERROR = 10e-20

def estimate_h(f, g, x):
    itr_count = 0
    h_estimated = [0]
    while True:
        cum = 0
        sigma_count = 0
        for sub_x in x:
            cum += ((g(sub_x) - f(sub_x + h_estimated[-1])) *
                    weight(sub_x) * f_1st_deriv(sub_x + h_estimated[-1]))
            sigma_count += weight(sub_x) * f_1st_deriv(sub_x + h_estimated[-1]) * f_1st_deriv(sub_x + h_estimated[-1])
        h_estimated.append(h_estimated[-1] + cum / sigma_count)

        if np.absolute(g(sub_x) - f(sub_x + h_estimated[-1])) < ERROR:
            break

        itr_count += 1
        if 100 == itr_count:
            break