2D Newton (in 2 Dimensionen) unter C# zur Findung von Lösungen und Nullstellen von Gleichungssystemen

2D Newton (in 2 Dimensionen) unter C# zur Findung von Lösungen und Nullstellen von Gleichungssystemen

Kleine (vor allem mathematische) Codestücke sind meistens nicht vorhanden, wenn man sie benötigt.

Alternativ sind sie in einer riesigen Bibliothek, wo die einzelnen Funktionen so stark verzahnt sind, daß man sie nicht auseinander reißen kann, und Dutzende MB extra einbinden muß – für eine Funktion.

Wenn man nur einen 2D Newton benötigt, hilft evtl. auch das hier.

Newton

1
2
3
4
5
6
7
8
9
10
11
12
13
       public static double[] Newton2Step(double[] x, Func<double, double, double> f, Func<double, double, double> g)
        {
            if (x.Length != 2)
                throw new ArgumentException("Should be 2x vector!");
            double[] res = new double[2];
            double[] temp = { f(x[0], x[1]), g(x[0], x[1]) };
            double[,] jacob = Jacobian2x2(x, f, g);
            Invert2x2Matrix(ref jacob);
            temp = Mult2x(jacob, temp);
            res[0] = x[0] - temp[0];
            res[1] = x[1] - temp[1];
            return res;
        }

Matrix Inverse

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
        static void Invert2x2Matrix(ref double[,] matrix)
        {
            // 0, 0 top left, 1, 1: top right
            if ((matrix.GetLength(0) != 2) || (matrix.GetLength(1) != 2))
            {
                throw new ArgumentException("Should be 2x2 matrix!");
            }
            double a, b, c, d;
            a = matrix[0, 0];
            b = matrix[1, 0];
            c = matrix[0, 1];
            d = matrix[1, 1];
            double factor = (1 / (a * d - b * c));
            matrix[0, 0] = factor * d;
            matrix[1, 0] = -factor * b;
            matrix[0, 1] = -factor * c;
            matrix[1, 1] = factor * a;
        }

Ableitung

1
2
3
4
5
6
7
8
9
        static double Derivative(double x, Func<double, double> f)
        {
            double delta = Math.Abs((x + 1E-3) * 1E-5);
            double x1 = x - delta;
            double x2 = x + delta;
            double y1 = f(x1);
            double y2 = f(x2);
            return (y2 - y1) / (x2 - x1);
        }

Jacobian Matrix

1
2
3
4
5
6
7
8
9
        static double[,] Jacobian2x2(double[] x, Func<double, double, double> f, Func<double, double, double> g)
        {
            double[,] res = new double[2, 2];
            res[0, 0] = Derivative(x[0], x_ => f(x_, 1));
            res[0, 1] = Derivative(x[1], x_ => f(1, x_));
            res[1, 0] = Derivative(x[0], x_ => g(x_, 1));
            res[1, 1] = Derivative(x[1], x_ => g(1, x_));
            return res;
        }

Matrizenmultiplikation und Vektormultiplikation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
        static double[,] Mult2x2(double[,] a, double[,] b)
        {
            double[,] res = new double[2, 2];
            res[0, 0] = a[0, 0] * b[0, 0] + a[0, 1] * b[1, 0];
            res[0, 1] = a[0, 0] * b[0, 1] + a[0, 1] * b[1, 1];
            res[1, 0] = a[1, 0] * b[0, 0] + a[1, 1] * b[1, 0];
            res[1, 1] = a[1, 0] * b[0, 1] + a[1, 1] * b[1, 1];
            return res;
        }
 
        static double[] Mult2x(double[,] a, double[] v)
        {
            double[] res = new double[2];
            res[0] = a[0, 0] * v[0] + a[0, 1] * v[1];
            res[1] = a[0, 0] * v[0] + a[0, 1] * v[1];
            return res;
        }

Sag es durch Code. 🙂 Dem ist kaum noch etwas hinzuzufügen, ich setze es hier mal als Referenz für alle Suchenden rein.

Referenzen

Das alles funktioniert gemäß der beim MIT gefundenen Infos.
Eine schöne Powerpoint zu dem Thema ist leider schon wieder offline.

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert.