関数の速度（4） - inamori’s diary

log(1+x)/xを有理関数で近似するとき、テーラー展開で4次まであわせると、

$\frac{\log{(1+x)}}{x} \approx \frac{30+21x+x^2}{30+36x+9x^2}$
$\qquad\qquad\qquad\qquad\qquad = 1-\frac{1}{2}x+\frac{1}{3}x^2-\frac{1}{4}x^3+\frac{1}{5}x^4-\frac{33}{200}x^5 + \cdots$

となったが、4次まであわせるには、

$\frac{1+a_1x}{1+b_1x+b_2x^2+b_3x^3}$

という形などが考えられる。
その中から最適なものを選ぶべきではないだろうか。

一般に、近似したい関数を、

$1 + c_1x + c_2x^2 + \cdots + c_nx^n + \cdots$

として、n次までの近似を求めたいとき、

$\frac{1 + a_1x + \cdots + a_lx^l}{1 + b_1x + \cdots + b_mx^m}$
ただし、l + m = n

が考えられる。
最初の例は連分数から有理関数を得たが、この考え方で求めてみよう。

$1 + a_1x + \cdots = (1 + b_1x + \cdots + b_mx^m)(1 + c_1x + \cdots)$
$a_k = c_k + b_1c_{k-1} + b_2c_{k-2}$
ただし、 $c_0 \equiv 1$

このとき、 $a_3 = a_4 = 0$ でなければならないから、

$c_3 + b_1c_2 + b_2c_1 = 0$
$c_4 + b_1c_3 + b_2c_2 = 0$
$c_k = \frac{(-1)^k}{k+1}$

（↑これ思いつくのに数日かかった）
この連立方程式、

$\frac{1}{3}b_1 - \frac{1}{2}b_2 = \frac{1}{4}$
$-\frac{1}{4}b_1 + \frac{1}{3}b_2 = -\frac{1}{5}$

を解いて、

$b_1 = \frac{6}{5} \qquad b_2 = \frac{3}{10}$
$a_1 = c_1 + b_1c_0 = \frac{7}{10}$
$a_2 = c_2 + b_1c_1 + b_2c_0 = \frac{1}{30}$

となり、上の結果と一致した。
一般には、

$a_{l+1} = c_{l+1}b_0 + \cdots + c_{l+1-m}b_m$
...
$a_{n} = c_{n}b_0 + \cdots + c_{n+1-m}b_m$
ただし、c_k = 0 (k < 0)

これは、連立一次方程式だから、b₁, ..., b_mについて解ける。
さらに、

$a_1 = c_1b_0 + \cdots + c_{1-m}b_m$
...
$a_l = c_lb_0 + \cdots + c_{l-m}b_m$

より、a₁, ..., a_lについても求められる。
求められた多項式をh_lと書くと、

$h_0(x) = \frac{1}{1+\frac{1}{2}x-\frac{1}{12}x^2+\frac{1}{24}x^3-\frac{19}{720}x^4}$
$h_1(x) = \frac{1+\frac{19}{30}x}{1+\frac{17}{15}x+\frac{7}{30}x^2-\frac{1}{90}x^3}$
$h_2(x) = \frac{1+\frac{7}{10}x+\frac{1}{30}x^2}{1+\frac{6}{5}x+\frac{3}{10}x^2}$
$h_3(x) = \frac{1+\frac{3}{10}x-\frac{1}{15}x^2+\frac{1}{60}x^3}{1+\frac{4}{5}x}$
$h_4(x) = 1-\frac{1}{2}x+\frac{1}{3}x^2-\frac{1}{4}x^3+\frac{1}{5}x^4$

ソースは一部のみ。



#include 

#include "taylor.h"
typedef vector    Vec;

typedef vector         Matrix;

typedef Taylor    Poly;
void solve(int n, int l);

fraction getTerm(int i);

void make_matrix(Matrix& mat, int n, int m);

void solve_matrix(Matrix& mat, Vec& b, int n);

void calc_a(Vec& a, const Vec& b, int l, int m);
void main(int argc, char **argv) {

    if(argc != 2) {

        cerr << "usage : log4 n." << endl;

        exit(1);

    }

    int     n = atoi(argv[1]);

    

    for(int l = 0; l <= n; l++)

        solve(n, l);

}
void solve(int n, int l) {

    int     m = n - l;

    Vec     a(l+1);

    Vec     b(m+1);

    Matrix  mat(m);

    make_matrix(mat, n, m);

    

    solve_matrix(mat, b, n);

    calc_a(a, b, l, m);

    a[0] = b[0] = 1;

    

    Poly    x("x");

    Poly    f, g;

    for(int i = l; i >= 0; i--)

        f = f * x + a[i];

    for(int i = m; i >= 0; i--)

        g = g * x + b[i];

    

    cout << f << endl;

    cout << g << endl;

//  cout << f / g << endl;

}
fraction getTerm(int i) {

    return fraction(i & 1 ? -1 : 1, i + 1);

}
void make_matrix(Matrix& mat, int n, int m) {

    int     l = n - m;

    for(int i = 0; i < m; i++) {

        mat[i].resize(m + 1);

        for(int k = 0; k < m; k++) {

            int     index = i - k + l;

            mat[i][k] = index >= 0 ? getTerm(index) : 0;

        }

        mat[i][m] = getTerm(l + i + 1) * -1;

    }

}
void solve_matrix(Matrix& mat, Vec& b, int n) {

    int     m = (int)b.size() - 1;

    

    // 上三角に

    for(int i = 0; i < m; i++) {

        if(mat[i][i].getBunsi() == 0) {

            // 対角が0にならないものを探してスワップ

            for(int j = i + 1; j < m; j++) {

                if(mat[j][i].getBunsi() != 0) {

                    Vec     tmp = mat[i];

                    mat[i] = mat[j];

                    mat[j] = tmp;

                    break;

                }

            }

        }

        

        // 対角を1に

        for(int k = i + 1; k <= m; k++)

            mat[i][k] /= mat[i][i];

        mat[i][i] = 1;

        

        // 対角成分より下を0に

        for(int j = i + 1; j < m; j++) {

            for(int k = i + 1; k <= m; k++)

                mat[j][k] -= mat[j][i] * mat[i][k];

        }

    }

    

    // 残りの非対角成分も0に

    for(int i = 1; i < m; i++) {

        for(int j = 0; j < i; j++) {

            for(int k = i + 1; k <= m; k++)

                mat[j][k] -= mat[i][k] * mat[j][i];

        }

    }

    

    // 解をコピー

    for(int i = 0; i < m; i++)

        b[i+1] = mat[i][m];

}
void calc_a(Vec& a, const Vec& b, int l, int m) {

    for(int i = 1; i <= l; i++) {

        a[i] = getTerm(i);

        for(int k = 1; k <= m; k++) {

            int     index_c = i - k;

            if(index_c >= 0)

                a[i] += getTerm(index_c) * b[k];

        }

    }

}