微分は引き算

java 数値計算

微分は引き算ってことがわかりやすくなるようなプログラムを考えてみる。
微分というのは、関数の傾きを求める操作です。


で、これを式として求めるときの規則として簡単なものに次のようなものがあります。
(x^n)'\rightarrow nx^{n-1}
例えば次のような式を微分します。
f(x)=10x^3+3x^2-x


そうすると、こうなります。xの0乗は1なので、最後の項からはxがなくなってます。
f'(x)=30x^2+6x-1


微分の定義としては、つぎのようになります。
f'(x)=\lim_{\Delta\to\0}\frac{f(x+\Delta)-f(x)}{\Delta}


これは、Δを無限に0に近づけたときの\frac{f(x+\Delta)-f(x)}{\Delta}ということになるのですが、プログラムで計算するときは無限に0に近づけるということはできないので、ほどほどに小さい値を使います。こういった求め方を差分といいます。
ここでは、式として求めた理論的な値と、差分の値がどんだけ違うか比べてみます。


メソッドfが目的の関数、dfが理論的な微分式で、それと並べて\frac{f(x+dx)-f(x)}{dx}を表示しています。
見た目、理論値と差分の値はほとんど変わりません。
精度が欲しいときには\frac{f(x+dx)-f(x-dx)}{2dx}を使うこともあります。

import java.awt.*;
import java.awt.image.BufferedImage;
import javax.swing.*;

public class Differentiation {
    //目的の関数
    private static double f(double x){
        return 10 * x * x * x - 3 * x * x - x;
    }
    //微分の理論値
    private static double df(double x){
        return 30 * x * x - 6 * x  - 1;
    }
    
    public static void main(String[] args){
        Image img = new BufferedImage(400, 300, BufferedImage.TYPE_INT_RGB);
        Graphics2D g = (Graphics2D) img.getGraphics();
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, 400, 300);
        g.setColor(Color.BLACK);
        g.drawString("目的関数 f(x) = 10 x^3 - 3 x^2 - x", 10, 20);
        g.drawString("理論値   f'(x) = 30 x^2 - 6 x - 1", 10, 160);
        g.drawString("差分     (f(x + dx) - f(x)) / dx", 10, 260);
        double dx = .01;
        for(int i = 0; i < 300; ++i){
            double x = (i - 150) * dx;
            //目的関数
            g.fillRect(i + 10, (int)(50 -  f(x)),                      2, 2);
            //微分の理論値
            g.fillRect(i + 10, (int)(150 - df(x)),                     2, 2);
            //差分
            g.fillRect(i + 10, (int)(250 - ((f(x + dx) - f(x)) / dx)), 2, 2);
        }
        g.dispose();
        
        JFrame f = new JFrame("微分");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setSize(400, 350);
        f.add(new JLabel(new ImageIcon(img)));
        f.setVisible(true);
    }
}

差分の数学的意味と誤差

数値計算

数学では、なんとなく同じような値になるからといって、無条件に使えるような都合のいい話はゆるされません。
差分が微分の近似として使えるのにもちゃんとした理由があります。
ここで、テイラー展開を持ってきます
\sum_{n=0}^\infty\frac{f^{(n)}(a)}{n!}(x-a)^n


ここで、x=x0+d、a=x0と置くと、次のようになります。
f(x_0+d)=f(x_0)+df'(x_0)+\frac{d^2}{2}f''(x_0)+\frac{d^3}{6}f'''(x_0)+\cdots
この式を変形すると
f'(x_0)=\frac{f(x_0+d)-f(x_0)}{d}-(\frac{d}{2}f''(x_0)+\frac{d^2}{6}f'''(x_0)+\cdots)
となります。dを十分小さい数値とするので
O(d)=\frac{d}{2}f''(x_0)+\frac{d^2}{6}f'''(x_0)+\cdots
は、誤差として切り捨てると、前進差分の式
f'(x_0)\approx\frac{f(x_0+d)-f(x_0)}{d}
が残ります。
つまり、この場合は、O(d)が誤差になります。


また、x=x0-d, a=x0と置いてテイラー展開すると次のようになります。
f(x_0-d)=f(x_0)-df'(x_0)+\frac{d^2}{2}f''(x_0)-\frac{d^3}{6}f'''(x_0)+\cdots
これをx=x0+d, a=x0と置いたテイラー展開から引くと
f(x_0+d)-f(x_0-d)=2df'(x_0)+2\frac{d^3}{6}f'''(x_0)+\cdots
となって、これを変形すると
f'(x_0)=\frac{f(x_0+d)-f(x_0-d)}{2d}-(\frac{d^2}{6}f'''(x_0)+\cdots)
となり、
O(d^2)=\frac{d^2}{6}f'''(x_0)+\cdots
を誤差として切り捨てると、中心差分の式
f'(x_0)\approx\frac{f(x_0+d)-f(x_0-d)}{2d}
が残ります。
この場合は、O(d^2)が誤差になります。


このとき、dは十分小さい値で、d>>d^2となるので、O(d)>>O(d^2)となります。
このことから、中心差分の方が前進差分より誤差が小さいといえることがわかります。

ついでに、2階微分の差分式

数値計算

上のふたつのテイラー展開を足すと
f(x_0-d)+f(x_0+d)=2f(x_0)+d^2f''(x_0)+2\frac{d^4}{4!}f''''(x_0)+\cdots
となります。
※f'''が消えたのでいままで省略されてたf''''をちゃんと書きました。
これを変形すると
f''(x_0)=\frac{f(x_0-d)-2f(x_0)+f(x_0+d)}{d^2}-(\frac{d^2}{12}f''''(x_0)+\cdots)
となります。
O(d^2)=\frac{d^2}{12}f''''(x_0)+\cdots
を誤差として切り捨てると、2階微分の差分式
f''(x_0)\approx\frac{f(x_0-d)-2f(x_0)+f(x_0+d)}{d^2}
が出てきます。


追記：
これを使って、ラプラス方程式など２階微分方程式を解くことができます。
ラプラシアンの計算をしてラプラス方程式を解く