Zenn
😺

JavaScriptで数式微分を作る(1)

2022/12/23に公開
JavaScript
math
微分
自動微分
tech

はじめに

授業の課題で p5.js を使って作品を作るという課題が出たので、どうせならと思い数式微分を作ることにしました。
冬休み中に作成し、3学期に提出できるようにします。

数式微分とは

Google で「数式微分」と調べてもそれについて解説している記事が見つかりませんでした。
Wiki でそれっぽいのを英語で見つけたのでおいておきます。
https://en.wikipedia.org/wiki/Computer_algebra

数値微分、自動微分との違い

このふたつは調べるとたくさん出てきます。

数値微分

https://en.wikipedia.org/wiki/Numerical_differentiation

数値微分は微分係数の定義に沿って近似するものです。
微分係数の定義は次の式であることはご存じの方も多いのではないでしょうか。

f^\prime\left(x\right)=\lim_{h\rightarrow0}{\frac{f\left(x+h\right)-f\left(x\right)}{h}}

手計算では極限をとることはできますが、プログラムではとることはできません。代わりにhを十分に小さい数にすることで近似します。

自動微分

https://ja.wikipedia.org/wiki/自動微分

自動微分（じどうびぶん、アルゴリズム的微分とも）とは、プログラムで定義された関数を解析し、偏導関数の値を計算するプログラムを導出する技術である。自動微分は複雑なプログラムであっても加減乗除などの基本的な算術演算や基本的な関数（指数関数・対数関数・三角関数など）のような基本的な演算の組み合わせで構成されていることを利用し、これらの演算に対して連鎖律を繰り返し適用することによって実現される。自動微分を用いることで偏導関数値を少ない計算量で自動的に求めることができる。
https://ja.wikipedia.org/wiki/自動微分

自動微分は連鎖律を用いています。例えば、y=f\left(g\left(x\right)\right)=f\left(w\right)では次のようになります。

\frac{dy}{dx}=\frac{dy}{dw}\frac{dw}{dx}

これにより簡単な微分を繰り返して元の式を微分します。

数式微分は最後に数値を代入するのに対し、この二つはすべての計算に数値を代入しているという違いがあるのではないでしょうか。

構想

最終的に数式微分を目標としますが、まずは自動微分を作ってそれを基に数式微分を作りたいと思います。
ユーザーが式を代入するときには WolframAlpha の数学代入みたいなものを使って入力できるようにしたいです。
https://ja.wolframalpha.com/

まずは自動微分

それではまず自動微分を作っていきます。
こちらのサイトを参考にして製作しました。
https://kivantium.hateblo.jp/entry/2016/03/25/010320

BUAD.js
/**
 * @typedef {Decimal | number | string} numbers
 */
class BUAD {
    /** @type {Decimal} */
    #val;
    /** @type {Decimal} */
    #dval;
    /** @type {Decimal.Config} */
    #decimalConfig;
    /** @type {Decimal.Constructor} */
    #BUADDecimal = Decimal;
    /**
     *
     * @param {numbers} v
     * @param {numbers} dv
     * @param {Decimal.Config|undefined} decimalConfig
     */
    constructor(v = 0, dv = 0, decimalConfig) {
        this.setBUADDecimal(decimalConfig);
        this.val = v;
        this.dval = dv;
    }

    selectX() {
        this.dval = new this.BUADDecimal("1");
    }

    /** @returns {Decimal} */
    get val() { return this.#val; }
    /** @returns {Decimal} */
    get dval() { return this.#dval; }
    get decimalConfig() { return this.#decimalConfig; }
    get BUADDecimal() { return this.#BUADDecimal; }

    /**
     * @param {numbers} v
     */
    set val(v) {
        if (v.constructor === Number || v.constructor === String) {
            this.#val = new this.BUADDecimal(v);
        } else if (v instanceof Decimal) {
            this.#val = new this.BUADDecimal(v.toString());
        } else {
            throw TypeError('Invalid value type.');
        }
    }

    /**
     * @param {numbers} dv
     */
    set dval(dv) {
        if (dv.constructor === Number || dv.constructor === String) {
            this.#dval = new this.BUADDecimal(dv);
        } else if (dv instanceof Decimal) {
            this.#dval = new this.BUADDecimal(dv.toString());
        } else {
            throw TypeError('Invalid value type.');
        }
    }

    /**
     * @param {Decimal.Config} config
     */
    set decimalConfig(config) { this.setBUADDecimal(config); }

    /**
     *
     * @param {Decimal.Config|undefined} object
     */
    setBUADDecimal(object) {
        if (object === undefined) return;
        this.decimalConfig = object;
        this.#BUADDecimal = this.#BUADDecimal.clone(object);
    }

    /**
     *
     * @param {numbers} v
     * @param {numbers} dv
     * @returns
     */
    create(v, dv) {
        return new BUAD(v, dv, this.decimalConfig);
    }

    copy() {
        return this.create(new this.BUADDecimal(this.val), new this.BUADDecimal(this.dval));
    }

    toString() {
        return `BUAD(${this.val.toString()}, ${this.dval.toString()})`;
    }

    // Four arithmetic operations and others
    /**
     *
     * @param {BUAD|undefined} x
     * @returns {BUAD}
     */
    add(x = undefined) {
        if (x === undefined) return this.copy();
        return this.create(this.val.add(x.val), this.dval.add(x.dval));
    }
    /**
     *
     * @param {BUAD|undefined} x
     * @returns {BUAD}
     */
    sub(x = undefined) {
        if (x === undefined) return this.create(this.val.negated(), this.dval.negated());
        return this.create(this.val.sub(x.val), this.dval.sub(x.dval));
    }

    /**
     *
     * @returns {BUAD}
     */
    negate() {
        return this.sub();
    }

    /**
     *
     * @param {BUAD} x
     * @returns {BUAD}
     */
    mul(x) {
        return this.create(this.val.mul(x.val), this.dval.mul(x.val).add(this.val.mul(x.dval)));
    }

    /**
     *
     * @param {BUAD} x
     * @returns {BUAD}
     */
    div(x) {
        return this.create(this.val.div(x.val), (this.dval.mul(x.val).sub(this.val.mul(x.dval))).div(x.val.pow(2)));
    }

    // Comparison functions
    /**
     *
     * @param {BUAD} x
     * @returns {boolean}
     */
    lt(x) {
        return this.val.lt(x.val);
    }

    /**
     *
     * @param {BUAD} x
     * @returns {boolean}
     */
    leq(x) {
        return this.val.lte(x.val);
    }

    /**
     *
     * @param {BUAD} x
     * @returns {boolean}
     */
    gt(x) {
        return this.val.gt(x.val);
    }

    /**
     *
     * @param {BUAD} x
     * @returns {boolean}
     */
    geq(x) {
        return this.val.gte(x.val);
    }

    //basic functions
    /**
     *
     * @returns {BUAD}
     */
    sqrt() {
        const sqrtX = this.val.sqrt();
        return this.create(sqrtX, this.dval.div(new this.BUADDecimal("2").mul(sqrtX)));
    }

    /**
     * pow(x, y)
     * return x^y, x^(y-1)(yx'+xy' ln(x))
     * @param {BUAD} y
     * @returns {BUAD}
     */
    pow(y) {
        const xPym1 = this.val.pow(y.val.sub(new this.BUADDecimal(1)));
        const yxP = y.val.mul(this.dval);
        const xyPlnx = (this.val.comparedTo(new Decimal(0)) === 0 || y.dval.comparedTo(new Decimal(0)) === 0) ? new Decimal(0) : this.val.mul(y.dval).mul(this.val.ln());
        return this.create(this.val.pow(y.val), xPym1.mul(yxP.add(xyPlnx)));
    }

    /**
     *
     * @returns {BUAD}
     */
    exp() {
        const expX = this.val.exp();
        return this.create(expX, this.dval.mul(expX));
    }

    /**
     *
     * @returns {BUAD}
     */
    log() {
        return this.create(this.val.ln(), this.dval.div(this.val));
    }

    /**
     *
     * @returns {BUAD}
     */
    sin() {
        return this.create(this.val.sin(), this.dval.mul(this.val.cos()));
    }

    /**
     *
     * @returns {BUAD}
     */
    cos() {
        return this.create(this.val.cos(), this.dval.mul(this.val.sin().negated()));
    }

    /**
     *
     * @returns {BUAD}
     */
    asin() {
        return this.create(this.val.asin(), this.dval.div(new this.BUADDecimal(1).sub(this.val.pow(2))));
    }

    /**
     *
     * @returns {BUAD}
     */
    acos() {
        return this.create(this.val.acos(), this.dval.negated().div(new this.BUADDecimal(1).sub(this.val.pow(2))));
    }

    /**
     *
     * @returns {BUAD}
     */
    atan() {
        return this.create(this.val.atan(), this.dval.div(new this.BUADDecimal(1).add(this.val.pow(2))));
    }
}

解説

このプログラムは Decimal.js を使って計算しています。
https://github.com/MikeMcl/decimal.js/

/**
 * @typedef {Decimal | number | string} numbers
 */

JsDoc の型定義です。numbers は BUAD の初期化時の引数に使います。v,dv の 2 つの型となります。

constructor周り
/** @type {Decimal} */
#val;
/** @type {Decimal} */
#dval;
/** @type {Decimal.Config} */
#decimalConfig;
/** @type {Decimal.Constructor} */
#BUADDecimal = Decimal;

/**
 *
 * @param {numbers} v
 * @param {numbers} dv
 * @param {Decimal.Config|undefined} decimalConfig
 */
constructor(v = 0, dv = 0, decimalConfig) {
    this.setBUADDecimal(decimalConfig);
    this.val = v;
    this.dval = dv;
}

/** @returns {Decimal} */
get val() { return this.#val; }
/** @returns {Decimal} */
get dval() { return this.#dval; }
get decimalConfig() { return this.#decimalConfig; }
get BUADDecimal() { return this.#BUADDecimal; }

/**
 * @param {numbers} v
 */
set val(v) {
    if (v.constructor === Number || v.constructor === String) {
        this.#val = new this.BUADDecimal(v);
    } else if (v instanceof Decimal) {
        this.#val = new this.BUADDecimal(v.toString());
    } else {
        throw TypeError('Invalid value type.');
    }
}

/**
 * @param {numbers} dv
 */
set dval(dv) {
    if (dv.constructor === Number || dv.constructor === String) {
        this.#dval = new this.BUADDecimal(dv);
    } else if (dv instanceof Decimal) {
        this.#dval = new this.BUADDecimal(dv.toString());
    } else {
        throw TypeError('Invalid value type.');
    }
}

/**
 * @param {Decimal.Config} config
 */
set decimalConfig(config) { this.setBUADDecimal(config); }

/**
 *
 * @param {Decimal.Config|undefined} object
 */
setBUADDecimal(object) {
    if (object === undefined) return;
    this.decimalConfig = object;
    this.#BUADDecimal = this.#BUADDecimal.clone(object);
}

微分前の値 val,微分後の値 dvalです。また、Decimal をこの中でクローンしているので Config と型を保持しておきます。

/**
 *
 * @param {numbers} v
 * @param {numbers} dv
 * @returns
 */
create(v, dv) {
    return new BUAD(v, dv, this.decimalConfig);
}

copy() {
    return this.create(new this.BUADDecimal(this.val), new this.BUADDecimal(this.dval));
}

toString() {
    return `BUAD(${this.val.toString()}, ${this.dval.toString()})`;
}

少し必要なやつです。

selectX() {
    this.dval = new this.BUADDecimal("1");
}

変数とするものにこれを適用します。

四則演算

足し算

add
/**
 *
 * @param {BUAD|undefined} x
 * @returns {BUAD}
 */
add(x = undefined) {
    if (x === undefined) return this.copy();
    return this.create(this.val.add(x.val), this.dval.add(x.dval));
}

引数が空の場合は符号が変わらないのでコピーを返します。
引数がある場合は微分前と微分後でそれぞれの和を返します。

f\left(x\right)=g\left(x\right)+h\left(x\right)\\ f^\prime\left(x\right)=g^\prime\left(x\right)+h^\prime\left(x\right)

引き算

sub
/**
 *
 * @param {BUAD|undefined} x
 * @returns {BUAD}
 */
sub(x = undefined) {
    if (x === undefined) return this.create(this.val.negated(), this.dval.negated());
    return this.create(this.val.sub(x.val), this.dval.sub(x.dval));
}

引数が空の場合は符号を反転します。
引数がある場合は微分前と微分後でそれぞれの差を返します。

f\left(x\right)=g\left(x\right)-h\left(x\right)\\ f^\prime\left(x\right)=g^\prime\left(x\right)-h^\prime\left(x\right)
negate
/**
 *
 * @returns {BUAD}
 */
negate() {
    return this.sub();
}

符号を反転するこれもあります。これもsubを引数なしで呼んでいるだけです。

掛け算

multiply
/**
 *
 * @param {BUAD} x
 * @returns {BUAD}
 */
mul(x) {
    return this.create(this.val.mul(x.val), this.dval.mul(x.val).add(this.val.mul(x.dval)));
}

微分前は積を、微分後は積の微分を返します。

f\left(x\right)=g\left(x\right)h\left(x\right)\\ f^\prime\left(x\right)=g^\prime\left(x\right)h\left(x\right)+g\left(x\right)h^\prime\left(x\right)

割り算

divide
/**
 *
 * @param {BUAD} x
 * @returns {BUAD}
 */
div(x) {
    return this.create(this.val.div(x.val), (this.dval.mul(x.val).sub(this.val.mul(x.dval))).div(x.val.pow(2)));
}

微分前は商を、微分後は商の微分を返します。

f\left(x\right)=\frac{g\left(x\right)}{h\left(x\right)}\\ f^\prime\left(x\right)=\frac{g^\prime\left(x\right)h\left(x\right)-g\left(x\right)h^\prime\left(x\right)}{\left(h\left(x\right)\right)^2}

比較演算子

/**
 *
 * @param {BUAD} x
 * @returns {boolean}
 */
lt(x) {
    return this.val.lt(x.val);
}

/**
 *
 * @param {BUAD} x
 * @returns {boolean}
 */
leq(x) {
    return this.val.lte(x.val);
}

/**
 *
 * @param {BUAD} x
 * @returns {boolean}
 */
gt(x) {
    return this.val.gt(x.val);
}

/**
 *
 * @param {BUAD} x
 * @returns {boolean}
 */
geq(x) {
    return this.val.gte(x.val);
}

それぞれ微分前の数値を比較しているだけです。

基本的な関数

平方根

squareRoot
/**
 *
 * @returns {BUAD}
 */
sqrt() {
    const sqrtX = this.val.sqrt();
    return this.create(sqrtX, this.dval.div(new this.BUADDecimal("2").mul(sqrtX)));
}
f\left(x\right)=\sqrt{g\left(x\right)}\\ f^\prime\left(x\right)=\frac{g^\prime\left(x\right)}{2\sqrt{g\left(x\right)}}

べき乗

power
/**
 * pow(x, y)
 * return x^y, x^(y-1)(yx'+xy' ln(x))
 * @param {BUAD} y
 * @returns {BUAD}
 */
pow(y) {
    const xPym1 = this.val.pow(y.val.sub(new this.BUADDecimal(1)));
    const yxP = y.val.mul(this.dval);
    const xyPlnx = (this.val.comparedTo(new Decimal(0)) === 0 || y.dval.comparedTo(new Decimal(0)) === 0) ? new Decimal(0) : this.val.mul(y.dval).mul(this.val.ln());
    return this.create(this.val.pow(y.val), xPym1.mul(yxP.add(xyPlnx)));
}
f\left(x\right)=g\left(x\right)^{h\left(x\right)}\\ f^\prime\left(x\right)=g\left(x\right)^{h\left(x\right)-1}\left(g^\prime\left(x\right)h\left(x\right)+g\left(x\right)h^\prime\left(x\right)\ln{g\left(x\right)}\right)

xyPlnxはg\left(x\right)h^\prime\left(x\right)\ln{g\left(x\right)}の部分を表しています。
g\left(x\right)=0の時にthis.val.ln()が NaN となって計算できなくなります。期待されるxyPlnxは 0 なので判定してあげてます。

exponential
/**
 *
 * @returns {BUAD}
 */
exp() {
    const expX = this.val.exp();
    return this.create(expX, this.dval.mul(expX));
}

これはf\left(x\right)=e^{g\left(x\right)}の時の計算です。

f\left(x\right)=e^{g\left(x\right)}\\ f^\prime\left(x\right)=g^\prime\left(x\right)e^{g\left(x\right)}

自然対数

logarithm
/**
 *
 * @returns {BUAD}
 */
log() {
    return this.create(this.val.ln(), this.dval.div(this.val));
}
f\left(x\right)=\ln{g\left(x\right)}\\ f^\prime\left(x\right)=\frac{g^\prime\left(x\right)}{g\left(x\right)}

三角関数

sin
/**
 *
 * @returns {BUAD}
 */
sin() {
    return this.create(this.val.sin(), this.dval.mul(this.val.cos()));
}
f\left(x\right)=\sin{g\left(x\right)}\\ f^\prime\left(x\right)=g^\prime\left(x\right)\cos{g\left(x\right)}
cos
/**
 *
 * @returns {BUAD}
 */
cos() {
    return this.create(this.val.cos(), this.dval.mul(this.val.sin().negated()));
}
f\left(x\right)=\cos{g\left(x\right)}\\ f^\prime\left(x\right)=g^\prime\left(x\right)\left(-\sin{g\left(x\right)}\right)

逆三角関数

arcsin
/**
 *
 * @returns {BUAD}
 */
asin() {
    return this.create(this.val.asin(), this.dval.div(new this.BUADDecimal(1).sub(this.val.pow(2))));
}
f\left(x\right)=\arcsin{g\left(x\right)}\\ f^\prime\left(x\right)=\frac{g^\prime\left(x\right)}{1-\left(g\left(x\right)\right)^2}
arccos
/**
 *
 * @returns {BUAD}
 */
acos() {
    return this.create(this.val.acos(), this.dval.negated().div(new this.BUADDecimal(1).sub(this.val.pow(2))));
}
f\left(x\right)=\arccos{g\left(x\right)}\\ f^\prime\left(x\right)=-\frac{g^\prime\left(x\right)}{1-\left(g\left(x\right)\right)^2}
arctan
/**
 *
 * @returns {BUAD}
 */
atan() {
    return this.create(this.val.atan(), this.dval.div(new this.BUADDecimal(1).add(this.val.pow(2))));
}
f\left(x\right)=\arctan{g\left(x\right)}\\ f^\prime\left(x\right)=\frac{g^\prime\left(x\right)}{1+\left(g\left(x\right)\right)^2}

実際に動かす

汚いですが次のコードで動かしました。

sketch.js
const canvasSize = 600;

function setup() {
  createCanvas(canvasSize, canvasSize);
  textAlign(RIGHT, TOP);
  textSize(20);
  drawFunction();
}
function drawFunction() {
  clear();
  background(255);
  grid();
  const x = new BUAD();
  x.val = 5;
  x.selectX();
  strokeWeight(2);
  for (let x = -500; x <= 500; x += 1) {
    const decimalX = new BUAD(Decimal(x.toString()).div(100));
    decimalX.selectX();
    // 半径5の円の上半分
    //const y=new BUAD(25).sub(decimalX.pow(new BUAD(2))).sqrt()
    // e^x^x
    //const y=decimalX.pow(decimalX).exp()
    // x^2*e^x
    //const y=decimalX.pow(new BUAD(2)).mul(decimalX.exp())
    // (e^x+e^-x)/2 cosh
    const y = decimalX.exp().add(decimalX.sub().exp()).div(new BUAD(2));
    // cos
    //const y = decimalX.cos();
    // x/2
    //const y = decimalX.div(new BUAD(2));
    // x^2
    //const y=decimalX.pow(new BUAD(2))
    //atan(asin(x))
    //const y = decimalX.asin().atan();

    stroke("black");
    point(
      ((5 + x / 100) * canvasSize) / 10,
      (-(-5 + y.val.toNumber()) * canvasSize) / 10
    );
    stroke("green");
    point(
      ((5 + x / 100) * canvasSize) / 10,
      (-(-5 + y.dval.toNumber()) * canvasSize) / 10
    );
  }
}
function grid() {
  strokeWeight(0.5);
  stroke("black");
  text("o", canvasSize / 2, canvasSize / 2);
  line(0, canvasSize / 4, canvasSize, canvasSize / 4);
  line(0, canvasSize / 2, canvasSize, canvasSize / 2);
  line(0, (canvasSize * 3) / 4, canvasSize, (canvasSize * 3) / 4);

  line(canvasSize / 4, 0, canvasSize / 4, canvasSize);
  line(canvasSize / 2, 0, canvasSize / 2, canvasSize);
  line((canvasSize * 3) / 4, 0, (canvasSize * 3) / 4, canvasSize);
}
y=\cosh\left(x\right)=\frac{e^x+e^{-x}}{2}

cosh

まとめ

今回はボトムアップ型自動微分を実装してみました。次回は数式微分に行く前にテキストで正しい式が出せるか試してみます。

Discussion