nichtlineares Gleichungssystem lösen

[Frank Buss]

Ich bin dabei, die Parameter für eine Kettenlinie auszurechnen. Ich bin
zu diesen Formeln gekommen:

h1 = a*cosh((0 - x0) / a) + d - a
h2 = a*cosh((w - x0) / a) + d - a

Die Unbekannten sind x0 und a, alle anderen Variablen sind bekannt.
Beispiel:

w = 10
h1 = 2.3769571985548783
h2 = 4.100221828554815
d = 2

Ergebnis:

a = 12
x0 = 3

Ich vermute mal, a und x0 kann man nicht symbolisch lösen? Bei
WolframAlpha zumindest kommt nur ein Timeout. Numerisch wäre aber auch
in Ordnung. Wie kann ich das machen?

Habe mal geschaut, geht wahrscheinlich mit dem Gauss-Newton Verfahren,
aber ich kenne nur Schulmathematik, keine Ahnung was eine Jacobi-Matrix
ist. Geht es einfacher?

Eine Beschreibung des Algorithmus wäre schön. Ich will das im Endeffekt
in C programmieren, wobei ich in Python teste, da ich mit matplotplib
recht schnell Tests zeichnen kann.

[Ralf Goertz]

Am Fri, 30 Apr 2021 12:08:07 +0200

Post by Frank Buss
Ich bin dabei, die Parameter für eine Kettenlinie auszurechnen. Ich
h1 = a*cosh((0 - x0) / a) + d - a
h2 = a*cosh((w - x0) / a) + d - a
Die Unbekannten sind x0 und a, alle anderen Variablen sind bekannt.
w = 10
h1 = 2.3769571985548783
h2 = 4.100221828554815
d = 2
a = 12
x0 = 3
Ich vermute mal, a und x0 kann man nicht symbolisch lösen? Bei
WolframAlpha zumindest kommt nur ein Timeout. Numerisch wäre aber
auch in Ordnung. Wie kann ich das machen?
Habe mal geschaut, geht wahrscheinlich mit dem Gauss-Newton
Verfahren, aber ich kenne nur Schulmathematik, keine Ahnung was eine
Jacobi-Matrix ist. Geht es einfacher?
Eine Beschreibung des Algorithmus wäre schön. Ich will das im
Endeffekt in C programmieren, wobei ich in Python teste, da ich mit
matplotplib recht schnell Tests zeichnen kann.

Ich würde ja davon abraten, das selber zu programmieren, denn trivial
ist das nicht. Die Jacobi-Matrix ist die Matrix der Ableitungen, bei
dir, da du 2 Variablen hast, ist es eine 2x2-Matrix. Aber man kommt ohne
aus. Ich würde dafür aber wie gesagt, etwas professionelles empfehlen.
Mit ein paar wenigen Zeilen geht es zum Beispiel in R
<https://www.r-project.org/>:

Definiere eine Funktion, die du minimieren möchtest. Du hast zwei
Gleichungen, die jeweils einen Wert ausspucken. Du hast Werte
vorgegeben, also schaust du dir die quadratische Abweichung der Werte
von den Funktionswerten mit gegebenen Parametern a und x0 an und
summierst die:

f=function(p,w=10,d=2,h1=2.3769571985548783,h2=4.100221828554815) {
a=p[1];
x0=p[2];
return( (a*cosh(-x0/a)+d-a-h1)^2 + (a*cosh((w-x0)/a)+d-a-h2)^2)
}

(Die beiden Parameter, die du berechnen willst, müssen in ein Feld hier
p, die restlichen Parameter habe ich mit in die Funktion gesteckt, aber
mit default-Werten belegt, weil sie ja bei jedem Aufruf gleich bleiben.)

Dann rufst du die Funktion nlm (non-linear minimisation) auf, deren
ersten beiden Parameter die zu minimierende Funktion und das Feld der zu
berechnenden Parameter sind. Dafür vergibst du zunächst Startwerte, ich
habe hier a=10 und x0=4 genommen. Die gibt dir dann im besten Fall das
Ergebnis (in estimate):

nlm(f,c(10,4))
$minimum
[1] 1.770767e-12

$estimate
[1] 11.999995 2.999999

$gradient
[1] -7.203196e-08 -2.882011e-07

$code
[1] 1

$iterations
[1] 12

Code 1 sagt dir, dass wohl alles glatt gelaufen ist. Und tatsächlich ist
das Ergebnis sehr nahe an (12; 3).

[Frank Buss]

Post by Ralf Goertz
f=function(p,w=10,d=2,h1=2.3769571985548783,h2=4.100221828554815) {
a=p[1];
x0=p[2];
return( (a*cosh(-x0/a)+d-a-h1)^2 + (a*cosh((w-x0)/a)+d-a-h2)^2)
}

R sieht interessant aus, werde ich mir bei Gelegenheit mal anschauen.
Kann ich aber im aktuellen Projekt nicht brauchen, da das C ist. Ich
habe aber die entsprechenden Funktionen aus der R-Library
zusammenkopiert, das meiste war hier zu finden:

https://github.com/wch/r-source/blob/trunk/src/appl/uncmin.c

Ein paar Fortran-Funktionen habe ich noch nach C übertragen. Ist jetzt
insgesamt eine 2500 Zeilen lange C-Datei, ohne Fortran Bindings. Damit
kann ich es jetzt auch berechnen. Also ganz einfach :-)

[Alfred Flaßhaar]

Post by Frank Buss
Ich bin dabei, die Parameter für eine Kettenlinie auszurechnen. Ich bin
h1 = a*cosh((0 - x0) / a) + d - a
h2 = a*cosh((w - x0) / a) + d - a
Die Unbekannten sind x0 und a, alle anderen Variablen sind bekannt.

(...)

Mathcad löst dieses kleine nichtlineare System wie ein Blitz so schnell.
Auch größere Abweichungen der Startwerte vom Ergebnis für die
Iterationen stören kaum.

Gruß, Alfred Flaßhaar

[Alfred Flaßhaar]

(...)

p. s.

Falls zur Veranschaulichung des math./phys. Hintergrundes besteht, ist
zu empfehlen:

Istvan Szabo, Technische Mechanik

Einführung, Band 1, ab S. 206
Repertorium ..., Band 3, ab S. 34

[Dieter Heidorn]

Post by Frank Buss
h1 = a*cosh((0 - x0) / a) + d - a
h2 = a*cosh((w - x0) / a) + d - a
Die Unbekannten sind x0 und a, alle anderen Variablen sind bekannt.
w = 10
h1 = 2.3769571985548783
h2 = 4.100221828554815
d = 2
a = 12
x0 = 3
[...] Numerisch wäre aber auch in Ordnung. Wie kann ich das machen?
Habe mal geschaut, geht wahrscheinlich mit dem Gauss-Newton Verfahren,
aber ich kenne nur Schulmathematik, keine Ahnung was eine Jacobi-
Matrix ist.

Du hast zwei Gleichungen, die von den zwei Variablen a und x0 abhängen.
Die schreibst du mit Funktionen in der Form

f1(a,x0) = 0 , f2(a,x0) = 0.

Die Jacobi-Matrix enthält die partiellen Ableitungen dieser Funktionen
nach den Variablen:

--- ---
| df1 df1 |
| --- --- |
| da dx0 |
J(a,x0) = | |
| df2 df2 |
| --- --- |
| da dx0 |
--- ---

Post by Frank Buss
Eine Beschreibung des Algorithmus wäre schön.

In Maxima sieht das Ganze so aus:

/* Konstanten: */;

w : 10 ; d : 2;
h1 : 2.3769571985548783 ; h2 : 4.100221828554815;

/* Variable: a, x0 */;

/* Funktionen: */;

f1(a,x0) := a*cosh((0 - x0) / a) + d - a - h1;
f2(a,x0) := a*cosh((w - x0) / a) + d - a - h2;

/* Jacobi-Matrix: */;

jacobian( [f1(a,x0),f2(a,x0)], [a,x0] );

--- ---
| x0*sinh(x0/a) |
| - ------------- + cosh(x0/a) - 1 sinh(x0/a) |
| a |
| |
| sinh((10-x0)/a)*(10-x0) sinh((10-x0 )|
| - ---------------------- + cosh((10-x0)/a) - 1 - ----------- |
| a a ---
---

load(mnewton);

/* Gleichungen: */

G1: f1(a,x0) = 0;
G2: f2(a,x0) = 0;

/* numerische Lösung: */;

/* Startwerte: a = 1, x0 = 1 */;

mnewton([G1,G2], [a,x0], [1,1]);

[[a=12.0,x0=3.0]]

Post by Frank Buss
Ich will das im Endeffekt in C programmieren,

Dann müsstest du das Newton-Verfahren beherrschen - in Maxima ist es in
der Funktion mnewton realisiert.

Wenn du zum Newton-Verfahren Fragen hast, stell' sie ruhig.

Dieter Heidorn

[Ernst Sauer]

Ich bin dabei, die Parameter für eine Kettenlinie auszurechnen. Ich bin zu diesen Formeln
h1 = a*cosh((0 - x0) / a) + d - a (1)
h2 = a*cosh((w - x0) / a) + d - a (2)
w = 10
h1 = 2.3769571985548783
h2 = 4.100221828554815
d = 2
a = 12
x0 = 3

Man kann auch die zweite Gleichung nach x0 lösen bzw. umformen
und dann in die erste Gleichung einsetzen.

Diese Gleichung löst man dann numerisch, z.B. nach Newton_Raphson.

Damit umgeht man die Jacobi-Matrix.

es

[Frank Buss]

Post by Ernst Sauer
Man kann auch die zweite Gleichung nach x0 lösen bzw. umformen

Und was kommt dabei raus? Ich hatte es mal mit WolframAlpha probiert,
und kam nichts Sinnvolles bei raus, hat nur noch eine weitere Unbekannte
erzeugt:

https://www.wolframalpha.com/input/?i=solve+h2+%3D+a*cosh%28%28w+-+x0%29+%2F+a%29+%2B+d+-+a%2Cx0

[Ernst Sauer]

Post by Ernst Sauer
Man kann auch die zweite Gleichung nach x0 lösen bzw. umformen

Und was kommt dabei raus? Ich hatte es mal mit WolframAlpha probiert, und kam nichts
https://www.wolframalpha.com/input/?i=solve+h2+%3D+a*cosh%28%28w+-+x0%29+%2F+a%29+%2B+d+-+a%2Cx0

Da steht ja auch nicht die umgeformte Gleichung.
Ermittle aus der zweiten Gleichung x0, also
x0 = ...

Diese Größe setzt Du dann in die erste Gleichung ein,
dort steht dann nur noch als Unbekannte die Größe a.

Für WolframAlpha benennst Du a um in x und gibst diese Gleichung
mit den bekannten Zahlenwerten für die anderen Größen ein.

Für die Nullstelle von y=a+b*x würdest du nur a+b*x eingeben

Bin aber kein Spezialist für WolframAlpha

Du wolltest doch ein C-Programm schreiben.
Für Newton_Raphson findest Du bestimmt den Code.

es

[Dieter Heidorn]

Post by Ernst Sauer

Post by Frank Buss
Ich bin dabei, die Parameter für eine Kettenlinie auszurechnen. Ich
h1 = a*cosh((0 - x0) / a) + d - a (1)
h2 = a*cosh((w - x0) / a) + d - a (2)
w = 10
h1 = 2.3769571985548783
h2 = 4.100221828554815
d = 2

Man kann auch die zweite Gleichung nach x0 lösen bzw. umformen
und dann in die erste Gleichung einsetzen.

Radio Eriwan meint: Im Prinzip ja... ;-)

Dabei ergibt sich aber ein Problem:

Formt man die zweite Gleichung schrittweise um, so kommt man zunächst
zu:

h2 - d + a ( w - x0 )
---------- = cosh(--------)
a ( a )

Der Hyperbelkosinus ist eine gerade Funktion - also müsste man zur
Umkehrung wissen, welches Vorzeichen das Argument des cosh hat.
Da man a und x0 nicht kennt, kann an dieser Stelle nichts dazu gesagt
werden.

Lässt man die Auflösung der zweiten Gleichung nach x0 von einem CAS
ausführen, so erhält man Meldungen wie z.B. in Maxima:

solve: using arc-trig functions to get a solution.
Some solutions will be lost.

Dieter Heidorn

[Ernst Sauer]

...

Post by Dieter Heidorn

Post by Ernst Sauer
Man kann auch die zweite Gleichung nach x0 lösen bzw. umformen
und dann in die erste Gleichung einsetzen.

Radio Eriwan meint: Im Prinzip ja... ;-)
Formt man die zweite Gleichung schrittweise um, so kommt man zunächst
h2 - d + a       ( w - x0 )
---------- = cosh(--------)
    a            (   a    )
Der Hyperbelkosinus ist eine gerade Funktion -  also müsste man zur
Umkehrung wissen, welches Vorzeichen das Argument des cosh hat.
Da man a und x0 nicht kennt, kann an dieser Stelle nichts dazu gesagt
werden.
Lässt man die Auflösung der zweiten Gleichung nach x0 von einem CAS
   solve: using arc-trig functions to get a solution.
   Some solutions will be lost.

Radio Eriwan meldet aber auch schon bei der Wurzel einer quadratischen
Gleichung, dass man aufpassen muss.

Der Anwender kennt sein Problem (hoffentlich) und muss beurteilen,
welche Lösung für das Problem brauchbar ist.
Er muss ja auch die Grenzen angeben in welchem das numerische Verfahren
die Nullstelle suchen soll und mit welcher Genauigkeit.
Das kann schon auch Fummelarbeit werden.
Die Kettenlinie ist keine simple Funktion, bzw. bei der Vorgabe
der Randbedingungen muss man schon aufpassen.

Ich würde mir auch immer den Graph der Lösungsfunktion anschauen.

es

[Alfred Flaßhaar]

...

Wozu diese Umstände? Viele CAS (z. B. Mathcad) haben starke numerische
Verfahren (Levenberg-Marquardt, konjugierte Gradienten usw.)
implementiert. Für diese gegenüber Startwerten robusten Verfahren ist
das numerische Lösen dieses kleinen Systems Sekundensache. Das Fahrrad
muß also nicht neu erfunden werden.

Freundliche Wochenendgrüße, Alfred Flaßhaar

[Ernst Sauer]

Post by Alfred Flaßhaar

...

Wozu diese Umstände?

Welche Umstände?

Post by Alfred Flaßhaar
Viele CAS (z. B. Mathcad) haben starke numerische Verfahren
(Levenberg-Marquardt, konjugierte Gradienten usw.) implementiert. Für diese gegenüber
Startwerten robusten Verfahren ist das numerische Lösen dieses kleinen Systems
Sekundensache. Das Fahrrad muß also nicht neu erfunden werden.

Die Lösung ist mit so einem Programm sicher eine Minutensache,
sofern man es hat (gekauft hat) und eingearbeitet ist.
Aber wenn das nicht gegeben ist, was dann?

Und da ist mein Vorschlag kurz und schmerzlos umsetzbar und für
den Anwender, der sein praktisches Problem im Griff hat und die
Schulmathematik beherrscht, auch sicher

Post by Alfred Flaßhaar
Freundliche Wochenendgrüße, Alfred Flaßhaar

Danke, wünsch ich ebenso.
E.S.