Lösen eines linearen Gleichungssystems (4x4)

[jakob78]

Hallo liebe Forianer,

ich stehe gerade vor einem Problem, dass ich mir einen Lösungsalgorithmus für ein lineares Gleichungssystem schreiben muss (in einer Fortran77-basierten Sprache). Und da ich durchs Mitlesen in diesem Forum hier schon mitbekommen habe, dass hier doch einige Programmierer/Mathematiker o.ä. anwesend sind, dachte ich, ich stelle mein Problem hier kurz einmal vor - eventuell hat der eine oder andere einen (oder mehrere) Hinweise für mich.

Nun handelt es sich dabei ein 4x4-Gleichungssystem, mit zum Teil sehr "ungeschickt" skalierten Matritzen, wie (wobei A = Koeffizientenmatrix und b = Störvektor) z.B.:

A = [ 3.91126953294906e-13  0.000290530407713016  0.118719995662381     0.999861140092325;
     -7.47034897265332e-13 -0.000158174186504138 -0.0169130360585405   -9.28246904902663e-06;
      1.42680307002891e-12  8.61151625166872e-05  0.00240945754016845   8.61761980650444e-11;
     -2.72512972030656e-12 -4.68838903437727e-05 -0.000343255085472549 -8.00038984641107e-16 ]

b = [ -0.000272415115940619; -9.83039417624391e-07;  1.42984225458907e-07;  3.29280620635692e-09 ]

Hierfür (aus Matlab: rank(A) = rank([A, b]) = 4, daher eindeutig lösbar) habe ich ein bestehendes Skript (Gauss-Elimination mit Pivoting) modifiziert bzw. in diese Fortran77-basierte Sprache portiert. Das funktioniert soweit auch ganz gut, wenngleich zunaechst die Matrix etwas skaliert werden muss, um die Zahlen dem float-Format etwas zugänglicher zu machen.

Das Problem tut sich nun auf, dass sich die Einträge des Gleichungssystem sukzessive verändern, bis irgendwann ein folgender Gestalt gelöst werden muss:

AA = [ 9.71932695068104e-21  1.95542531016539e-06  0.0319594622483047    0.999775498752868;
      -7.90937295921570e-21 -4.53931579859449e-07  -0.00194341016115096  -3.96425686254142e-06;
       6.43647249705794e-21  1.05375479248690e-07  0.000118176051434193  1.57188613761892e-11;
      -5.23785873026922e-21 -2.44618178588262e-08 -7.18612026002067e-06 -6.23275967050862e-17 ]

bb = [ 4.42616984297712e-05; -1.36475494887893e-06;  9.48249135038941e-08;  3.66966454607991e-09 ]

Hierfür erhalte ich in Matlab rank(AA) = 3 und rank([AA, bb]) = 4, das GLS ist also überbestimmt.

Nun stellen sich mir folgende Fragen:

- Ist das GLS AA,bb tatsächlich überbestimmt, oder ergibt sich dieses Problem aus der ungünstigen Skalierung der Matrix bzw. dieser schon "sehr kleinen" Einträge? (Eine zunächste Skalierung mit 10e+22 führt aber zum selben Ergebnis; weshalb das GLS tatsächlich überbestimmt sein dürfte.)

- Welches Verfahren eignet sich für die Lösung des GLS AA,bb. Ich dachte hier an eine QR-Zerlegung mit kleinsten Fehlerquadraten zum Finden des Lösungsvektors.

- Gibt es einen fertigen Sourcecode o.ä., mit dem beide Probleme behandelt/gelöst werden können? Ich bin leider kein Programmierer/Mathematiker und weiß daher nicht genau, ob ich selbst einen Löser so schreiben kann, dass dieser effizient und sicher die beiden GLS behandeln kann (noch dazu, wenn man auf das float-Zahlenformat beschränkt ist).

Über Rückmeldungen wäre ich sehr dankbar!

Liebe Grüße   jakob