C言語で偏微分のプログラムを実装

C言語で偏微分のプログラムを実装します。

定義

まずは偏微分の定義のおさらいから。

2変数関数$f(x,y)$は、$(a,b) \in \mathbb{R}^2$の周りで定義されているとする。このとき、

$$\begin{align*} f_x(a,b)&=\lim_{h \to 0}{f(a+h,b)-f(a,b) \over h} \end{align*}$$

が存在するならば，それを点$(a,b) \in \mathbb{R}^2$における$x$に関する偏微分 (partialderivative) という。
また、点$(a,b) \in \mathbb{R}^2$における$y$に関する偏微分

$$\begin{align*} f_x(a,b)&=\lim_{h \to 0}{f(a,b+h)-f(a,b) \over h} \end{align*}$$

も同様に定義される。

ソースコード

今回のコードでは数値微分の手法を用いることで偏微分を実現しています。
数値微分とは数値解析の分野の一つで，与えられる関数の常微分の値を解析的にではなく数値的に求める手法のことです。

PDef.c

#include <stdio.h>

/*任意の関数f(x,y)の定義*/
double f(double x, double y) 
{
    return 2*x*x + y*y + 5;
}

/*xに関する偏微分*/
double Fn_Def_x(double (*func)(double), double x, double y, double h) 
{
    return (func(x + h,y) - func(x,y)) / h;
}

/*yに関する偏微分*/
double Fn_Def_y(double (*func)(double), double x, double y, double h) 
{
    return (func(x, y + h) - func(x, y)) / h;
}

int main() 
{
    double h = 0.0001; // 微小な変化量   
    for (int x = 1; x < 4; x++)
    {
        for (int y = 1; y < 4; y++)
        {
            double def_x = Fn_Def_x(f, x, y, h);
            double def_y = Fn_Def_y(f, x, y, h);
            printf("∂f(%d,%d)/∂x＝ %f\n", x,y,def_x);
            printf("∂f(%d,%d)/∂y＝ %f\n", x,y,def_y);
        }
    }   
    return 0;
}

実行結果

実行結果は以下の通り。

出力結果

∂f(1,1)/∂x＝ 4.000200
∂f(1,1)/∂y＝ 2.000100
∂f(1,2)/∂x＝ 4.000200
∂f(1,2)/∂y＝ 4.000100
∂f(1,3)/∂x＝ 4.000200
∂f(1,3)/∂y＝ 6.000100
∂f(2,1)/∂x＝ 8.000200
∂f(2,1)/∂y＝ 2.000100
∂f(2,2)/∂x＝ 8.000200
∂f(2,2)/∂y＝ 4.000100
∂f(2,3)/∂x＝ 8.000200
∂f(2,3)/∂y＝ 6.000100
∂f(3,1)/∂x＝ 12.000200
∂f(3,1)/∂y＝ 2.000100
∂f(3,2)/∂x＝ 12.000200
∂f(3,2)/∂y＝ 4.000100
∂f(3,3)/∂x＝ 12.000200
∂f(3,3)/∂y＝ 6.000100

誤差はあるものの、どの点においても概ね正しく計算することができています。