C# で方程式を解く

// この投稿は C# Advent Calendar 2014 の 7 日目の記事です。

ゲーム アプリケーションやセンサーを用いたアプリケーションでは、
力学系の方程式の計算や座標系の変換など、何らかの数値計算が必要になることがあります。

中でも一次方程式は頻度が高いと思いますが、開発者の多くは毎回、
あらかじめ手計算により解の形式 (x = ~) を求めてから実装しているのではないでしょうか?
例えば、ax + b = cx + d の形式のものを x = (d – b) / (a – c) に変形してからプログラムするとか。

処理速度に重点を置く場合はこれでよいのですが、
最初に脳内で導出された式とソースコード上の式とでは形式に乖離が生じることがあり、
確認のために再計算しなければならなくなるなど、保守性が高いとはいえないでしょう。

パターン化できるのであれば、面倒な作業は自動化したいところです。
そこで以下では、なるべく元の方程式のままの形から解を得る仕組みを考えます。

 

まず、数式を C# の標準的な算術演算子を用いて表現できるようにするため、
一変数の多項式を表す Polynomial 構造体を以下のように定義します。
それぞれの次数に対する係数を Dictionary で管理し、必要な演算子を用意します。
(全体のソースコードは GitHub の EquationConsole にあります。)


public struct Polynomial
{
    public static readonly Polynomial X = new Polynomial(new Dictionary<int, double> { { 1, 1 } });

    static readonly IDictionary<int, double> _coefficients_empty = new Dictionary<int, double>();

    IDictionary<int, double> _coefficients;

    IDictionary<int, double> Coefficients
    {
        get { return _coefficients == null ? _coefficients_empty : _coefficients; }
    } 

    // The dictionary represents index/coefficient pairs.
    // The dictionary does not contain entries whose coefficient is 0.
    public Polynomial(IDictionary<int, double> coefficients)
    {
        _coefficients = coefficients;
    }

    public static implicit operator Polynomial(double value)
    {
        return value == 0 ? default(Polynomial) : new Polynomial(new Dictionary<int, double> { { 0, value } });
    }

    public static Polynomial operator +(Polynomial p1, Polynomial p2)
    {
        var coefficients = new Dictionary<int, double>(p1.Coefficients);

        foreach (var item2 in p2.Coefficients)
        {
            AddMonomial(coefficients, item2.Key, item2.Value);
        }
        return new Polynomial(coefficients);
    }

    public static Polynomial operator (Polynomial p1, Polynomial p2)
    {
        var coefficients = new Dictionary<int, double>(p1.Coefficients);

        foreach (var item2 in p2.Coefficients)
        {
            AddMonomial(coefficients, item2.Key, item2.Value);
        }
        return new Polynomial(coefficients);
    }

    public static Polynomial operator *(Polynomial p1, Polynomial p2)
    {
        var coefficients = new Dictionary<int, double>();

        foreach (var item1 in p1.Coefficients)
        {
            foreach (var item2 in p2.Coefficients)
            {
                AddMonomial(coefficients, item1.Key + item2.Key, item1.Value * item2.Value);
            }
        }
        return new Polynomial(coefficients);
    }

    public static Polynomial operator /(Polynomial p, double value)
    {
        var coefficients = new Dictionary<int, double>();

        foreach (var item in p.Coefficients)
        {
            AddMonomial(coefficients, item.Key, item.Value / value);
        }
        return new Polynomial(coefficients);
    }

    // Power
    public static Polynomial operator ^(Polynomial p, int power)
    {
        if (power < 0) throw new ArgumentOutOfRangeException("power", "The value must be non-negative.");

        Polynomial result = 1;

        for (var i = 0; i < power; i++)
        {
            result *= p;
        }
        return result;
    }

    public static Polynomial operator +(Polynomial p)
    {
        return p;
    }

    public static Polynomial operator (Polynomial p)
    {
        return 1 * p;
    }

    static void AddMonomial(Dictionary<int, double> coefficients, int index, double coefficient)
    {
        if (coefficients.ContainsKey(index))
        {
            coefficient += coefficients[index];
        }

        if (coefficient != 0)
        {
            coefficients[index] = coefficient;
        }
        else
        {
            coefficients.Remove(index);
        }
    }
}


今回は等値演算子 (==) などは省きました。
暗黙的型変換演算子 (implicit operator Polynomial) を定義したため、
多項式 -x + 2 や 2x – 1 を次のように表すことができます。

var x = Polynomial.X;
var p1 = x + 2;
var p2 = 2 * x  1;

 

次に、この Polynomial 構造体に、方程式を解くためのメソッドを追加します。
方程式の左辺はこの多項式自身で、右辺は 0 であるとします。
一次方程式のほか、二次方程式についても実装しています。


public int Degree
{
    get { return Coefficients.Count == 0 ? 0 : Coefficients.Max(c => c.Key); }
}

// Substitution
public double this[double value]
{
    get { return Coefficients.Sum(c => c.Value * Math.Pow(value, c.Key)); }
}

// Solve the equation whose right operand is 0.
public double SolveLinearEquation()
{
    if (Degree != 1) throw new InvalidOperationException("The degree must be 1.");

    // ax + b = 0
    var a = GetCoefficient(1);
    var b = GetCoefficient(0);

    return b / a;
}

// Solve the equation whose right operand is 0.
public double[] SolveQuadraticEquation()
{
    if (Degree != 2) throw new InvalidOperationException("The degree must be 2.");

    // ax^2 + bx + c = 0
    var a = GetCoefficient(2);
    var b = GetCoefficient(1);
    var c = GetCoefficient(0);
    var d = b * b  4 * a * c;

    return d > 0 ? new[] { (b  Math.Sqrt(d)) / (2 * a), (b + Math.Sqrt(d)) / (2 * a) }
        : d == 0 ? new[] { b / (2 * a) }
        : new double[0];
}

double GetCoefficient(int index)
{
    return Coefficients.ContainsKey(index) ? Coefficients[index] : 0;
}


インデクサーには代入の機能を割り当てました。
以上で Polynomial 構造体の最低限の実装はできたと思います。

これを使うと、例えば直線 y = x – 1 と直線 y = -2x + 5 の交点 P は次のように求められます。
ほぼ数式のままの形で記述できていると思います。


static class Program
{
    static readonly Polynomial x = Polynomial.X;

    static void Main(string[] args)
    {
        var l1 = x  1;
        var l2 = 2 * x + 5;
        var p_x = (l1 l2).SolveLinearEquation();
        var p_y = l1[p_x];
        Console.WriteLine("P({0}, {1})", p_x, p_y); // P(2, 1)
    }
}


 

次の例では、特定の 2 点を通る直線上で、y 座標から x 座標を求めます。


static class Program
{
    static readonly Polynomial x = Polynomial.X;

    static void Main(string[] args)
    {
        var p1 = new Point2D(0, 300);
        var p2 = new Point2D(1800, 300);
        var y_to_x = GetFunc_y_to_x(p1, p2);
        Console.WriteLine(y_to_x(0)); // 900
        Console.WriteLine(y_to_x(100)); // 600
    }

    // P1, P2 を通る直線上で、指定された y 座標に対応する x 座標を求めるための関数。
    static Func<double, double> GetFunc_y_to_x(Point2D p1, Point2D p2)
    {
        // P1(x1, y1) および P2(x2, y2) を通る直線の方程式 (の公式):
        // (x – x1) (y2 – y1) – (x2 – x1) (y – y1) = 0
        return y => ((x p1.X) * (p2.Y p1.Y) (p2.X p1.X) * (y p1.Y)).SolveLinearEquation();
    }
}

struct Point2D
{
    public double X { get; private set; }
    public double Y { get; private set; }

    public Point2D(double x, double y)
        : this()
    {
        X = x;
        Y = y;
    }
}


座標系のマッピングも、簡単なものであればこれとほぼ同様にできます (一般的には行列を使いますが)。

ちなみに、二次方程式 x2 + x -1 = 0 の解もこの通り。

二次方程式

 

なお、今回は引数が (Polynomial, Polynomial) のみの +, -, * 演算子を定義しましたが、
(Polynomial, double) などの引数を受け付けるオーバーロードを用意すれば、演算量を若干削減できると思います。

作成したサンプル
EquationConsole (GitHub)

参照
オーバーロードされた演算子 (C# プログラミング ガイド)

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中

%d人のブロガーが「いいね」をつけました。