Rückwärtsmultiplikation einer invertierten matrix

[Eichelhäer]

Hallo,

hab grade die inversen von 4x4 Matrizen berechnet (eigene Klasse usw. war viel schreibaufwand) und komme auch auf die richtigen werte verglichen mit diversen online Matrizen rechnern. Jetzt versuche ich den test zu machen gemäß der formel A*A^-1 = E und wider erwarten ist das Ergebnis unbefriedigend. Ich verwende double als werte und erhalte teilweise die form der einheitsmatrix. Neben den diagonalen einsen habe ich statt der übrigen nullen längere dezimalwerte und nicht die einheits matrix.
Hab schon versucht die nachkommastellen auf eine zu reduzieren : Math.round(10.0*(1/det))/10.0, aber erhalte nur nullen. Größere zehnerpotenzen ergeben wieder die Dezimalwerte in langer form. Hab auch schon versucht mit integern zu arbeiten hat aber nix gebracht, den bei 4x4 ist die det "meistens" kleiner 0. Wie bekomme ich die einheitsmatrix heraus in ihrer schönen diagonalform?

 

[httpdigest]

Da Fliesskommaarithmetik immer ungenau ist bzw. Rundungsungenauigkeiten hat, kannst du niemals sicherstellen, dass du den exakten mathematischen Wert herausbekommst. Du wirst hier also einfach mit Ungenauigkeiten leben müssen. Eventuell hilft dir java.lang.BigDecimal, aber auch hier wird bei nicht exakt als Dezimalzahl darstellbaren Brüchen (die du beim Invertieren bekommst) gerundet. Du kannst nur selber festlegen, ab welcher Stelle gerundet werden soll.
Eine weitere Möglichkeit wäre, deine Matrixelemente durch rationale Zahlen zu repräsentieren, also durch echte Brüche p/q und alle Operationen auf Brüchen zu implementieren.

 

[httpdigest]

Hier ist eine Implementierung mit rationalen Zahlen und BigInteger für Zähler und Nenner:

Spoiler

import java.math.BigInteger;
class Rational {
  BigInteger p, q;
  Rational(BigInteger p, BigInteger q) {
    BigInteger gdc = gdc(p, q);
    this.p = p.divide(gdc);
    this.q = q.divide(gdc);
  }
  private static BigInteger gdc(BigInteger a, BigInteger b) {
    while (!b.equals(BigInteger.ZERO)) {
      BigInteger tmp = b;
      b = a.mod(b);
      a = tmp;
    }
    return a;
  }
  static Rational d(double d) {
    long j = 1, num;
    do {
      j = j * 10;
    } while ((d * j) % 10 != 0);
    j = j / 10;
    num = (long) (d * j);
    return new Rational(BigInteger.valueOf(num), BigInteger.valueOf(j));
  }
  static Rational l(long p) {
    return new Rational(BigInteger.valueOf(p), BigInteger.ONE);
  }
  Rational mul(Rational b) {
    return new Rational(p.multiply(b.p), q.multiply(b.q));
  }
  Rational sub(Rational b) {
    return q == b.q ? new Rational(p.subtract(b.p), q) : new Rational(p.multiply(b.q).subtract(b.p.multiply(q)), q.multiply(b.q));
  }
  Rational add(Rational b) {
    return q == b.q ? new Rational(p.add(b.p), q) : new Rational(p.multiply(b.q).add(b.p.multiply(q)), q.multiply(b.q));
  }
  Rational negate() {
    return new Rational(p.negate(), q);
  }
  Rational reciprocal() {
    return new Rational(q, p);
  }
  public String toString() {
    if (q.equals(BigInteger.ONE))
      return p.toString();
    return p + "/" + q;
  }
}
public class RMatrix4x4 {
  Rational m00, m01, m02, m03;
  Rational m10, m11, m12, m13;
  Rational m20, m21, m22, m23;
  Rational m30, m31, m32, m33;
  public RMatrix4x4() {
    this(1, 0, 0, 0,
         0, 1, 0, 0,
         0, 0, 1, 0,
         0, 0, 0, 1);
  }
  public RMatrix4x4(
      long m00, long m01, long m02, long m03,
      long m10, long m11, long m12, long m13,
      long m20, long m21, long m22, long m23,
      long m30, long m31, long m32, long m33) {
    this(Rational.l(m00), Rational.l(m01), Rational.l(m02), Rational.l(m03),
         Rational.l(m10), Rational.l(m11), Rational.l(m12), Rational.l(m13),
         Rational.l(m20), Rational.l(m21), Rational.l(m22), Rational.l(m23),
         Rational.l(m30), Rational.l(m31), Rational.l(m32), Rational.l(m33));
  }
  public RMatrix4x4(
      double m00, double m01, double m02, double m03,
      double m10, double m11, double m12, double m13,
      double m20, double m21, double m22, double m23,
      double m30, double m31, double m32, double m33) {
    this(Rational.d(m00), Rational.d(m01), Rational.d(m02), Rational.d(m03),
         Rational.d(m10), Rational.d(m11), Rational.d(m12), Rational.d(m13),
         Rational.d(m20), Rational.d(m21), Rational.d(m22), Rational.d(m23),
         Rational.d(m30), Rational.d(m31), Rational.d(m32), Rational.d(m33));
  }
  public RMatrix4x4(
      Rational m00, Rational m01, Rational m02, Rational m03,
      Rational m10, Rational m11, Rational m12, Rational m13,
      Rational m20, Rational m21, Rational m22, Rational m23,
      Rational m30, Rational m31, Rational m32, Rational m33) {
    this.m00 = m00; this.m01 = m01; this.m02 = m02; this.m03 = m03;
    this.m10 = m10; this.m11 = m11; this.m12 = m12; this.m13 = m13;
    this.m20 = m20; this.m21 = m21; this.m22 = m22; this.m23 = m23;
    this.m30 = m30; this.m31 = m31; this.m32 = m32; this.m33 = m33;
  }
  public RMatrix4x4 invert() {
    Rational a = m00.mul(m11).sub(m01.mul(m10)), b = m00.mul(m12).sub(m02.mul(m10));
    Rational c = m00.mul(m13).sub(m03.mul(m10)), d = m01.mul(m12).sub(m02.mul(m11));
    Rational e = m01.mul(m13).sub(m03.mul(m11)), f = m02.mul(m13).sub(m03.mul(m12));
    Rational g = m20.mul(m31).sub(m21.mul(m30)), h = m20.mul(m32).sub(m22.mul(m30));
    Rational i = m20.mul(m33).sub(m23.mul(m30)), j = m21.mul(m32).sub(m22.mul(m31));
    Rational k = m21.mul(m33).sub(m23.mul(m31)), l = m22.mul(m33).sub(m23.mul(m32));
    Rational det = a.mul(l).sub(b.mul(k)).add(c.mul(j)).add(d.mul(i)).sub(e.mul(h)).add(f.mul(g)).reciprocal();
    return new RMatrix4x4(
      m11.mul(l).sub(m12.mul(k)).add(m13.mul(j)).mul(det),
      m01.negate().mul(l).add(m02.mul(k)).sub(m03.mul(j)).mul(det),
      m31.mul(f).sub(m32.mul(e)).add(m33.mul(d)).mul(det),
      m21.negate().mul(f).add(m22.mul(e)).sub(m23.mul(d)).mul(det),
      m10.negate().mul(l).add(m12.mul(i)).sub(m13.mul(h)).mul(det),
      m00.mul(l).sub(m02.mul(i)).add(m03.mul(h)).mul(det),
      m30.negate().mul(f).add(m32.mul(c)).sub(m33.mul(b)).mul(det),
      m20.mul(f).sub(m22.mul(c)).add(m23.mul(b)).mul(det),
      m10.mul(k).sub(m11.mul(i)).add(m13.mul(g)).mul(det),
      m00.negate().mul(k).add(m01.mul(i)).sub(m03.mul(g)).mul(det),
      m30.mul(e).sub(m31.mul(c)).add(m33.mul(a)).mul(det),
      m20.negate().mul(e).add(m21.mul(c)).sub(m23.mul(a)).mul(det),
      m10.negate().mul(j).add(m11.mul(h)).sub(m12.mul(g)).mul(det),
      m00.mul(j).sub(m01.mul(h)).add(m02.mul(g)).mul(det),
      m30.negate().mul(d).add(m31.mul(b)).sub(m32.mul(a)).mul(det),
      m20.mul(d).sub(m21.mul(b)).add(m22.mul(a)).mul(det)
    );
  }
  private static Rational mad(Rational a, Rational b, Rational c) {
    return a.mul(b).add(c);
  }
  public RMatrix4x4 mul(RMatrix4x4 m) {
    return new RMatrix4x4(
        mad(m00, m.m00, mad(m10, m.m01, mad(m20, m.m02, m30.mul(m.m03)))),
        mad(m01, m.m00, mad(m11, m.m01, mad(m21, m.m02, m31.mul(m.m03)))),
        mad(m02, m.m00, mad(m12, m.m01, mad(m22, m.m02, m32.mul(m.m03)))),
        mad(m03, m.m00, mad(m13, m.m01, mad(m23, m.m02, m33.mul(m.m03)))),
        mad(m00, m.m10, mad(m10, m.m11, mad(m20, m.m12, m30.mul(m.m13)))),
        mad(m01, m.m10, mad(m11, m.m11, mad(m21, m.m12, m31.mul(m.m13)))),
        mad(m02, m.m10, mad(m12, m.m11, mad(m22, m.m12, m32.mul(m.m13)))),
        mad(m03, m.m10, mad(m13, m.m11, mad(m23, m.m12, m33.mul(m.m13)))),
        mad(m00, m.m20, mad(m10, m.m21, mad(m20, m.m22, m30.mul(m.m23)))),
        mad(m01, m.m20, mad(m11, m.m21, mad(m21, m.m22, m31.mul(m.m23)))),
        mad(m02, m.m20, mad(m12, m.m21, mad(m22, m.m22, m32.mul(m.m23)))),
        mad(m03, m.m20, mad(m13, m.m21, mad(m23, m.m22, m33.mul(m.m23)))),
        mad(m00, m.m30, mad(m10, m.m31, mad(m20, m.m32, m30.mul(m.m33)))),
        mad(m01, m.m30, mad(m11, m.m31, mad(m21, m.m32, m31.mul(m.m33)))),
        mad(m02, m.m30, mad(m12, m.m31, mad(m22, m.m32, m32.mul(m.m33)))),
        mad(m03, m.m30, mad(m13, m.m31, mad(m23, m.m32, m33.mul(m.m33)))));
  }
  @Override
  public String toString() {
    return m00 + " " + m01 + " " + m02 + " " + m03 + "\n" +
           m10 + " " + m11 + " " + m12 + " " + m13 + "\n" +
           m20 + " " + m21 + " " + m22 + " " + m23 + "\n" +
           m30 + " " + m31 + " " + m32 + " " + m33 + "\n";
  }

  // Test
  public static void main(String[] args) {
    RMatrix4x4 m = new RMatrix4x4(
        0.2, 0, 0.63173, 5,
        0, 0.1, 0, 0,
        0, 0, 200, 2412023,
        2.4, -2.1643, 2, 1);
    System.out.println(m);
    RMatrix4x4 inv = m.invert();
    System.out.println(inv);
    System.out.println(m.mul(inv));
  }
}