オブジェクト指向で数値計算

はじめに

この記事は明治大学Advent Calender 2018の3日目の記事です。3日目にして2回目の投稿です。自転車操業Advent Calenderです。よろしくお願いします。

目標

• オブジェクト指向で数値計算を実装する。
• 数値計算を通してオブジェクト指向を理解する。

(注意1) 使用言語はProcessing(ほぼJava)である。Javaとの違いは描画面が充実している点である。
(注意2) オブジェクト指向で実装するから計算速度が早くなる、とかそんなことはない。メリット・デメリットを書いてある節もあるのでそちらを参照。

本記事の構成

• 数値計算法の紹介(Euler法):
  • Euler法を知っている方は飛ばしてok
• オブジェクト指向の用語:
  • クラスとインターフェースを知っている方は飛ばしてok
• オブジェクト指向で数値計算:
  • メインパート。読んでください。
• メリット・デメリット:
  • メインパートを読んだあとに読んでください。
• 発展的な話題:
  • オブジェクト指向の詳しい話を書きました。

数値計算法の紹介(Euler法)

Euler(オイラー)法の紹介をする。Euler法とはスキーム(数値計算の手法)の一つである。まずは一般的な形で紹介した後に、具体的な問題だとどうなるか紹介する。次のような常微分方程式の初期値問題が与えられたとする。

\frac{dx}{dt} = f(x, t), \quad x(0 ) = x_0

$f$の形によって手計算により解けるならそうすれば良いが、解けない場合は数値計算に頼るのが手っ取り早い。Euler法は上の常備分方程式を離散化して、以下のような漸化式を作って、解くことに相当する。なお、$dt$は時間の差分幅であり、どれだけ細かく数値計算をするかのパラメータとなる。

x_{n + 1} = x_n + dt \cdot f(x_n, \ dt \cdot n)

この漸化式によって逐次的に$x_n$を求めていけばよいのである。

簡単な例ではMalthus(マルサス)モデルという人口増加を表す常微分方程式がある。

\frac{dx}{dt} = ax, \quad x(0) = x_0

これは手計算でも$x(t) = x_0 e^{at}$という風に簡単に解くことができるが、Euler法により数値計算をするときは初期値を$x_0$として次のような漸化式を考えれば良い。

x_{n + 1} = x_n + dt \cdot (ax_n)

あとはこれを解けば、「何となく」本当の解に近そうな答えを得られる。

オブジェクト指向の用語

簡単にオブジェクト指向をおさらいしておく。オブジェクト指向とはざっくり言うと、「もの」や「概念」を一つの「クラス」として見るような考え方である。「クラス」において必要な事項を実装するように要請するものを「インターフェース」と呼ぶ。

よくある例を使って説明する。オブジェクト指向で「犬」を実装したいとする。我々が「犬」に要請したいことは例えば吠えることである。なので、「犬インターフェース」ではbark(吠える)関数を持つ。

ただ、一口に吠えると言っても犬種によって吠え方は違う。なので、「犬インターフェース」を継承した「プードルクラス」を作って実際に吠え方を「実装」してあげよう。

あとは、ペットとして飼いたいから、個々のプードルに名前をつけてあげたい。そうなったら、名前をもたせて上げればよいのである。個々のプードルのことを「犬クラスのインスタンス」と呼ぶ。

このアイディアをJavaで実装するとこうなる。

いぬ.java

interface Dog {
    // 継承するクラスに実装を要請する
    void bark();
}

class Poodle implements Dog {
    String name;

    Poodle (String name) {
        this.name = name;
    }

    // インターフェースの関数は全部実装する
    void bark () {
        System.out.println("きゃんきゃん");
    }
}

// Dog型としてPoodleのインスタンスを作ってる。
Dog alto = new Poodle("アルト");
alto.bark(); // きゃんきゃん

オブジェクト指向で数値計算

さて、本題である「オブジェクト指向で数値計算」であるが、MalthusモデルのEuler法を例にとってどう考えていくか見ていこう。先にイメージとしてはこんな感じである。本節を読んでいてわからなくなったら適宜下図に戻ってくると良い。

全体像

数値計算全体は大きく分けて「計算パート」と「可視化パート」に分かれる。「計算パート」は更に「モデルの設定」と「スキームの設定」がある。例えば今回ならMalthusモデルというモデルをEuler法で解く、といった具合である。以下で細かいところを見ていって、最終的に全体としてまとめ上げるようなボトムアップ方式で説明を行っていく。

Malthusモデル in 常微分方程式

Malthusモデルは無数にある常微分方程式の初期値問題の中の一つである。これをオブジェクト指向の言葉に言い換えれば、Malthusモデルというクラスは、常微分方程式の初期値問題というインターフェースを実装したものである。

なので、常微分方程式の初期値問題という大きなインターフェースとして「ODEインターフェース」を作ってあげよう。常微分方程式は$$\frac{dx}{dt} = f(x, t), \quad x(0) = x_0$$という風に一般的にかけるため、「ODEインターフェース」を実装するそれぞれのモデルに対して「$f$の実装」を要請する。

そして具体的なMalthusモデルを「ODEインターフェース」を実装することで作成するのである。Malthusモデルは人口増加率を表すパラメータ$a$を用いて、$f(x, t) = ax$と定義されるため、パラメータ$a$と初期値が変数となる。

以上を加味して実装すると次のようになる。

public interface ODE {
    public float getInitialVal ();
    public float f (float x, float t);
}


public class Malthus implements ODE {
    // パラメータ
    private float a;
    private float initialVal;

    public Malthus (float a, float initialVal) {
        this.a = a;
        this.initialVal = initialVal;
    }

    // ODEインターフェースの要請を実装
    public float getInitialVal () {
        return this.initialVal;
    }

    public float f (float x, float t) {
        return this.a * x;
    }
}

Euler法 in 数値計算スキーム

Euler法は沢山ある数値計算スキームの中の一つである。例えば常微分方程式の数値計算スキームは他にはRunge-Kutta法などがある。今回はそんな中からEuler法を選んだのである。これをオブジェクト指向の言葉で言うならば、Euler法というクラスは常微分方程式の数値計算スキームというインターフェースを実装したものである。

常微分方程式の数値計算スキームのインターフェースとして「Schemeインターフェース」を作成してあげよう。ここで各スキームが共通して持っていてほしい関数を考える。常微分方程式の数値計算スキームは漸化式を構成することだということを考えれば、次の2つの要請が考えられる。

• calcNextX(ODEモデル, 今の$x$, 今の$t$)
  • 与えられたモデルと、現在の$x$と$t$から次の時間$t + dt$での$x$の値を返す関数
• getSolution(ODEモデル)
  • calcNextXを逐次的に使用することで近似解をすべての離散時間で求める関数

あとは、「Schemeインターフェース」を実装した「Euler法クラス」を作成すれば良い。$dt$を時間の差分幅として、calcNextXは$$x_{n + 1} = x_n + dt \cdot f(x_n, n \cdot dt)$$で表される漸化式を実装すれば良い。

以上のことを実際にプログラムにしてみる。

public interface Scheme {
    public float calcNextX(ODE sys, float currentX, float currentT);
    public float[] getSolution (ODE sys);
}

public class Euler implements Scheme {
    private float dt;
    private int timeStep;

    public Euler (float dt, int timeStep) {
        this.dt = dt;
        this.timeStep = timeStep;
    }

    public float calcNextX(ODE sys, float currentX, float currentT) {
        return currentX + sys.f(currentX, currentT) * this.dt;
    }

    public float[] getSolution(ODE sys) {
        float[] xs = new float[this.timeStep];
        xs[0] = sys.getInitialVal();

        for (int i = 1; i < this.timeStep ; i++) {
             xs[i] = calcNextX(sys, xs[i - 1], (float)i * this.dt);
        }

        return xs;
    }
}

ODEとSchemeを統合

これまで、「ODEインターフェース」と「Schemeインターフェース」を独立に作ってきたことに気がついただろうか。「Scheme」では引数に「ODE」を取るような関数がありはしたが、「ODE」がどう実装されているかは一切気にしていない。

「ODE」と「Scheme」を一括して扱うにはもうひとつ上のレイヤーとして、「数値計算クラス」という大きな枠組みを作る必要がある。いわば管理者のような役割を持つクラスである。

内容は至極単純で、「ODE」と「Scheme」を一つずつ持って、あとは与えられた「ODE」と「Scheme」を用いて得た近似解を持っていればよいのである。

実装は以下のようになる。

public class NumericalAnalysis {
  private ODE sys;
  private Scheme scheme;
  private float[] solution;

  public NumericalAnalysis (ODE sys, Scheme scheme) {
    this.sys = sys;
    this.scheme = scheme;
    this.solution = this.getResult();
  }

  private float[] getResult() {
    return this.scheme.getSolution(this.sys);
  }
}

可視化

どうせなら数値計算結果を可視化したい。ここからはProcessing独特の書き方なので、詳細な説明は省くが、オブジェクト指向的な部分だけ説明する。

可視化も一つの概念であるため、「Graphicクラス」を作成する。これは数値計算にて得られた結果を解釈して描画を行うクラスである。

ここで注意するのが、あくまでも「Graphicクラス」の立場からすると、与えられたfloat配列を描画するだけであって、決して数値計算結果を描画するわけではない。オブジェクト指向では数値計算部分と可視化部分は完全に分業されるべきなのである。

一応かなりお粗末な可視化クラスをいかに載せておく。軸を書いて近似解の軌跡を描くだけである。もしグラフィックライブラリが使える環境なら是非をそっちを使ったほうが良い。

public class Graphic {
    public void drawAxis () {
       background(255);
       stroke(0);
       // x-axis
       line(30, 30, 30, 470);
       // t-axis
       line(30, 470, 470, 470);
    }

    public void drawGraph (float[] solution) {
        drawAxis();

        float maxVal = max(solution);
        float minVal = min(solution);
        float dt = 440 / (float)solution.length;

        for (int i = 0; i < solution.length - 1; i++) {
            float now = map(solution[i], minVal, maxVal, 465, 35);
            float next = map(solution[i + 1], minVal, maxVal, 465, 35);
            line(dt * i + 30, now, dt * (i + 1) + 30, next);
        }
    }
}

数値計算と可視化の統合

「数値計算クラス」は与えられた「ODE」と「Scheme」を用いて、数値計算結果を保存するようなクラスであった。また、「可視化クラス」はfloatの配列を描画するようなものであった。これらを統合して、「数値計算クラス」により得られた数値計算結果を「可視化クラス」に受け渡すようなクラスを作成する必要がある。

この役目を担うのが「Overallクラス」である。中身は単純で「数値計算クラス」と「可視化クラス」を結びつけて、数値計算から可視化までを一つの関数(doAll関数)を呼ぶだけでできるようにする。

public class Overall {
    private NumericalAnalysis na;
    private Graphic g;

    public Overall (ODE sys, Scheme scheme) {
        this.na = new NumericalAnalysis (sys, scheme);
        this.g = new Graphic();
    }

    public void doAll () {
        g.drawGraph(this.na.getResult());
    }
}

実行

以上であらかた完成したので、あとは実行のみ。MalthusモデルのインスタンスとEuler法のインスタンスをOverallに投げて、doAllにより実行すれば数値計算結果が可視化される。

void setup()はProcessing特有の書き方だが、main関数だと思ってもらえれば良い。

void setup() {
    // 画面サイズ
    size(500, 500);

    ODE sys = new Malthus(1.0, 1.0);
    Scheme scheme = new Euler(0.01, 300);
    Overall oa = new Overall (sys, scheme);

    oa.doAll();
}

実行結果

ちゃんと指数関数っぽい結果が得られていそう。(軸に目盛やラベルがないのは許してください。)

メリット・デメリット

オブジェクト指向で数値計算、どうだっただろうか。オブジェクト指向を使わなければもっと簡単に書けるのにという主張、その通りである。抽象化に抽象化を重ねたプログラムの構成をしたので、ちょっと複雑に見えるのかもしれない。

ただ、こうしてODEやスキームを抽象化したインターフェースを作成した結果、MalthusモデルやEuler法ではない別のモデルやスキームを作成することがた容易くなっているはずである。例えば、Euler法でなく4次のRunge-Kutta法で実装したいとなったら、「Schemeインターフェース」から「Runge-Kuttaクラス」を実装してそのインスタンスを作るだけで、もう可視化までできてしまうのである。つまり、1つの部品を交換するだけで別のプログラムが書けてしまうのである。

発展的な話題

ODE以外(PDEなど)

今回は常微分方程式の初期値問題とその数値計算スキームに的を絞って説明してきたが、他にも偏微分方程式とか、それに付随する色んな数値計算手法がある。ただ、かなりの枠組みが今回と同様の枠組みで説明できるため、一般化はさほど難しい話ではない。

SchemeをSingletonに

本記事ではスキームを普通のクラスとして書いた。これはすなわち、一個のプログラムの中でいくつも同じスキームのインスタンスを作れてしまうことを意味する。それはまずい。スキームはあくまでも方法であるため、同じものが複数あってはならないのだ。

それを解消する方法として、Singletonというアイディアがある。簡単に説明すると、プログラム内でそのクラスのインスタンスが高々1つのみしか生成されないことを保証するためのアイディアである。詳しいことはデザインパターン(リンクはWikipedia)を参照。

Singletonにすべき、というのは1例であり、リンク先のデザインパターンを見てもらえば、この記事のプログラムにどんなデザインパターンが使われているかとか、どれが使われるべきかとかがわかると思う。

MVC(Model・View・Controller)

「数値計算クラス」と「可視化クラス」の2つに分けたのは何でなのかと思った人はいるだろうか。可視化節で、「可視化クラス」はあくまでもfloat配列を翻訳して描画するだけで、数値計算結果を描画しているわけではないと強調して書いた。

この言葉の裏にはMVCという概念が潜んでいる。MVCとはModel View Controllerの略である。この概念はゲームを例に取るとわかりやすく、Modelがゲームのハードやソフト、Viewがゲームの画面、Controllerがコントローラーである。これらは人間や信号を介することで互いに繋がっている。画面に映る情報をもとに、コントローラーに入力し、その入力された情報をハードが解釈して、その結果を画面へと受け渡す。これらは全て互いに同じレベルで情報を共有しているのではなく、情報を翻訳して授受している。受け渡しのサイクルになってはいるが、それらを繋いでいるのは人間やケーブルなどの外部である。

今回の数値計算プログラムではControllerこそなかったが、数値計算クラスがModelであり、可視化クラスがViewであった。そのため、可視化クラスはfloat配列を描画するだけ、と言ったのである。

ポリモーフィズム

ポリモーフィズムとは日本語で多態性と訳される。今回、インターフェースとしてSchemeをインターフェースとして用意して、それをEuler法やRunge-Kutta法などが実装される枠組みを用意した。その中で、Euler法インスタンスを次のように作成した。

Scheme euler = new Euler(0.01, 1.0);

こうすることで、Scheme型なのにも関わらず、Euler法の動作を持つような変数を作ることができた。このように、抽象的な型のもとで実態としては具体的な型としていられるようなアイディアをポリモーフィズムと呼ぶ。ポリモーフィズムはプログラムの場合分けを1つにまとめられる役割などがある。

カプセル化

カプセル化とは端的に言えばブラックボックス化である。外から見れば中の煩雑な実装方法を気にすることなく、そのクラスを使うことができるようにする概念である。例えば今回は全てをまとめ上げたOverallクラスを用いてdoAll関数1つだけで全てが実行された。ユーザーはMalthusモデルとEuler法の各種パラメータを定めて、doAllするだけで中身を知ることなく結果が得られるのである。

まとめ

オブジェクト指向を用いて常微分方程式の数値計算を行うプログラムを作成した。そのプログラムを用いて、オブジェクト指向とはどういうものなのかを理解した。