PDA

Archiv verlassen und diese Seite im Standarddesign anzeigen : Längenbestimmung Hermitespline über Integration


Neomi
2007-03-25, 17:29:08
Eigentlich bin ich recht fit in Mathematik, aber momentan habe ich einen für mich unlösbaren Hirnknoten. Ich möchte die Länge eines beliebigen Hermitesplines exakt bestimmen können, ohne eine Unmenge von Iterationen zu benötigen, bei denen sich dann auch noch Rundungsfehler aufsummieren.

Konkret geht es um diese Funktion, für die ich das Integral (über u) nicht selbst bestimmen kann:

sqrt(A*u^4 + B*u^3 + C*u^2 + D*u + E)

Für das eigentliche Problem ist das zwar nicht relevant, aber die Parameter A bis E (alles Skalare) kommen so zustande...

A = va.x*va.x + va.y*va.y + va.z*va.z
B = 2*(va.x*vb.x + va.y*vb.y + va.z*vb.z)
C = 2*(va.x*vc.x + va.y*vc.y + va.z*vc.z) + vb.x*vb.x + vb.y*vb.y + vb.z*vb.z
D = 2*(vb.x*vc.x + vb.y*vc.y + vb.z*vc.z)
E = vc.x*vc.x + vc.y*vc.y + vc.z*vc.z

bzw.

A = dot(va,va)
B = 2*dot(va,vb)
C = 2*dot(va,vc) + dot(vb,vb)
D = 2*dot(vb,vc)
E = dot(vc,vc)

mit...

va = 3*(2*(p0-p1) + t0 + t1)
vb = 2*(3*(p1-p0) - 2*t0 - t1)
vc = t0

Die drei Vektoren va, vb und vc sind die Parameter, mit denen die Tangente des Hermitesplines (über p0, t0, p1, t1 definiert) bestimmt wird (t(u)=va*u^2 + vb*u + vc). Das oben angegebene Polynomial 4. Grades ist das Skalarprodukt der Tangente mit sich selbst, die Wurzel daraus ist die Länge der Tangente.

Über den Wolfram Integrator (http://integrals.wolfram.com/index.jsp) habe ich es schon versucht ("Sqrt[a*x^4 + b*x^3 + c*x^2 + d*x + f]", mit f statt e, da e als Eulersche Zahl vordefiniert ist), aber leider ohne Erfolg. Das Ergebnis ist kein "geht nüscht", sondern ein unfertiges Formelgedöhns mit Platzhaltern, das scheinbar vor der Zusammenkürzung zu lang und deshalb abgebrochen wurde. Google hat auch nichts ausgespuckt.

Gibt es zu dem Integral eine Lösung? Kann es überhaupt eine Lösung geben? Oder muß ich weiterhin iterativ rechnen?

Spasstiger
2007-03-25, 19:41:31
Mal ein Vorschlag: Mach eine Taylorreihenentwicklung von der Funktion bis die Entwicklung die gewünschte Genauigkeit erreicht und integriere dann die einzelnen Summanden der Reihe.
Für den Entwicklungspunkt 1 und eine Reihe bis einschließlich der vierten Ordnung komme ich mit Maple auf folgendes Polynom als Lösung des Integrals:

1/1920*(125*D^4-720*E^3*B+960*A^3*C*u+4200*B^2*u*C*D+11200*A^2*u^3*C*E+23520*A*B*E*D*u+27120*B*C *E*D*u-6792*A^2*B*D+24*B^3*C*u^5+2880*A*E^2*D*u^2+4800*A*E^2*C*u^4-6080*A*E^2*C*u^3+5760*A*E^2*C*u^2-168*A*C*D^2*u^5+2100*A*C*D^2*u^4-3360*A*C*D^2*u^3+10680*A*C*D^2*u^2-96*A*B*C^2*u^5+840*A*B*C^2*u^4+1120*A*B*C^2*u^3+8400*A*B*C^2*u^2+24*A*C*B^2*u^5-60*A*C*B^2*u^4+4000*A*C*B^2*u^3+5400*A*C*B^2*u^2+240*A*E*D^2*u^5+600*A*E*D^2*u^4-480*A*E*D^2*u^3+4080*A*E*D^2*u^2+192*C^2*E*D*u^5-1320*C^2*E*D*u^4-768*A*E^2*C*u^5-240*A^2*C^2*u-240*A*B^3*u+6720*A*E^3*u+8400*C^2*E^2*u+2160*C^2*D^2*u+4800*C^3*E*u+1800*C*D^3*u +2310*B^2*D^2*u+8400*B^3*E*u+4800*A^3*D*u-320*E^3*C+960*C^3*u*D+20160*A^2*C*u*E+3360*A*C^2*u*D+16800*C*E^2*u*D+1120*B^2*u^ 3*C^2+2800*B^2*u^3*A^2+4200*B^2*u^2*C^2+1680*B^2*u^2*A^2-240*B^2*u*A^2+2240*A^3*u^3*C-120*C*D*E^2*u^4+1680*C*D*E^2*u^2+240*A*C^3*u^4-60*C^2*B^2*u^4+192*A^2*B*D*u^5-15*D^4*u^5+105*D^4*u^4-350*D^4*u^3+1050*D^4*u^2+525*D^4*u+9*B^4*u^5-75*B^4*u^4+450*B^4*u^3+450*B^4*u^2-75*B^4*u+240*B*E^3*u^4-2016*A*B*E*D*u^5+10080*A*B*E*D*u^4-8960*A*B*E*D*u^3+9600*A*B*E*D*u^2+4000*C^2*E*D*u^3+1200*C^2*E*D*u^2+72*C*D^2*E*u ^5-420*C*D^2*E*u^4+1120*C*D^2*E*u^3+4200*C*D^2*E*u^2+840*B^2*C*E*u^5-5460*B^2*C*E*u^4+16800*B^2*C*E*u^3-14520*B^2*C*E*u^2-336*B^2*E*D*u^5+1050*B^2*E*D*u^4+2520*B^2*E*D*u^3-420*B^2*E*D*u^2+672*B*C^2*E*u^5-4200*B*C^2*E*u^4+12320*B*C^2*E*u^3-8400*B*C^2*E*u^2+8400*D^2*u*E^2+6720*D*u*E^3-960*A^2*E*u^4*C-2080*A^2*C^2-1494*A*B^3-640*A^4+2920*B^2*u^3*C*D+5460*B^2*u^2*C*D+1440*A^2*C^2*u^4-800*A^2*C^2*u^3+4800*A^2*C^2*u^2+24*A*B^3*u^5-150*A*B^3*u^4+1800*A*B^3*u^3+1500*A*B^3*u^2-1008*A^2*D^2*u^5+5400*A^2*D^2*u^4-720*A*C^3-1752*A^2*D^2-1792*A^2*E^2-1152*A*E^3-490*A*D^3-9720*A*B^2*E-210*C*B^2*u^4*D-832*C^2*E^2+188*C^2*D^2-752*C^3*E-1400*A*D^3*u^3+630*D^3*u^4*A+30*D^3*u^4*E+3360*u^2*B*C^3-240*A^2*C^2*u^5-1320*A^2*B*D*u^4+9600*A^2*B*D*u^3-5520*A^2*B*D*u^2+1680*A*B^2*E*u^5-10200*A*B^2*E*u^4+30240*A*B^2*E*u^3-31920*A*B^2*E*u^2-768*D*A^2*C*u^5+4200*D*A^2*C*u^4-3360*D*A^2*C*u^3+8400*D*A^2*C*u^2+1920*B*A^2*E*u^5-11760*B*A^2*E*u^4+35840*B*A^2*E*u^3-40320*B*A^2*E*u^2-96*A^2*B*C*u^5+600*A^2*B*C*u^4+4000*A^2*B*C*u^3+4080*A^2*B*C*u^2+240*A*B^2*D*u^5-320*B*C^3-279*B^4+1680*C*D*E*u^4*A+3360*C*D*E*u^2*A+16320*B*C^2*E*u+7200*B^2*E^2*u+40320*A *B*C*E*u+6720*C*E^3*u+4560*B*C*D^2*u+2100*B^3*D*u+1120*C*D*E^2*u^3-4632*D*A^2*C-12080*B*A^2*E+12960*B*D^2*E*u+14400*A^3*E*u-5224*A^2*B*C+298*C*D^3-5150*A*B^2*D-2850*A*B*D^2-4992*A^2*E*D-8320*A^2*C*E+8400*A*B*C*D*u+2240*A^3*u^3*B+640*A^4*u^3+320*B*C^3*u^3+320*C*E^3*u ^3-2360*A*C^2*D-4032*A*B*E^2-5152*A*C^2*E-3152*A*E^2*D-3712*A*E^2*C-1932*A*C*D^2-3064*A*B*C^2-336*A*B*C*D*u^5+2280*A*B*C*D*u^4+2240*A*B*C*D*u^3+11760*A*B*C*D*u^2-3964*A*C*B^2+16800*B*D*E^2*u+1860*B*D^3*u-2760*A*E*D^2+1120*B*D*u^3*C^2-9184*A*B*E*D-760*C*D*E^2-7280*A*C*D*E-1192*C^2*E*D-13440*A^2*u^2*C*E+8160*u^2*B*C^2*D-292*C*D^2*E-7064*A*B*C*D+480*A^3*u^2*B-80*D^2*u^3*E^2-280*D^3*u^3*E+1680*D^2*u^2*E^2-13904*A*B*C*E-5460*B^2*C*E-2290*B^2*C*D-3654*B^2*E*D+3360*B*u*C^2*D+16320*B*C*E^2*u-4904*B*C*E*D+80*E^2*D^2+190*E*D^3-8400*A^2*D^2*u^3+10080*A^2*D^2*u^2-2688*A^2*E^2*u^5+13440*A^2*E^2*u^4-22400*A^2*E^2*u^3+19200*A^2*E^2*u^2+192*A*E^3*u^5-48*C^2*E^2*u^5+1120*C^2*E^2*u^3-48*C^2*D^2*u^5+300*C^2*D^2*u^4-800*C^2*D^2*u^3+5400*C^2*D^2*u^2+192*C^3*E*u^5-1200*C^3*E*u^4+3200*C^3*E*u^3-1440*C^3*E*u^2-48*C*D^3*u^5+330*C*D^3*u^4-1000*C*D^3*u^3+3900*C*D^3*u^2-210*B^3*C*u^4+1400*B^3*C*u^3-3752*B*C^2*E-1120*B*C^2*D-742*B*C*D^2-2760*B*C*E^2+480*D*u^2*E^3+2100*B^3*C*u^2-600*B^3*D*u^4+2520*B^3*D*u^3-720*B^2*E^2*u^5+3240*B^2*E^2*u^4-3600*B^2*E^2*u^3+2880*B^2*E^2*u^2-210*B^2*D^2*u^5+1050*B^2*D^2*u^4-980*B^2*D^2*u^3+4500*B^2*D^2*u^2+504*B^3*E*u^5-3150*B^3*E*u^4+9000*B^3*E*u^3-8820*B^3*E*u^2-60*B*D^3*u^5+480*B*D^3*u^4-1000*B*D^3*u^3+3600*B*D^3*u^2+960*A^3*E*u^5-5760*A^3*E*u^4+16640*A^3*E*u^3-19200*A^3*E*u^2+192*A^3*D*u^5-1200*A^3*D*u^4+5440*A^3*D*u^3-4800*A^3*D*u^2+84*B^3*D*u^5+2240*A*D*u^3*C*E+8400*u^2*C^2*D*A+840*A*C^2*u^4*D-914*B^3*C-1224*B^3*D-1800*B^2*E^2-910*B^2*D^2-940*B^2*C^2-2574*B^3*E-1514*B*D^2*E-1944*B*D*E^2+1344*A*B*C*E*u^5-8400*A*B*C*E*u^4+31040*A*B*C*E*u^3-30240*A*B*C*E*u^2+144*B*C*E*D*u^5-1320*B*C*E*D*u^4+8000*B*C*E*D*u^3-240*B*C*E*D*u^2-168*B*C*D^2*u^5+1050*B*C*D^2*u^4-1400*B*C*D^2*u^3+9660*B*C*D^2*u^2-480*B*C*E^2*u^5+2040*B*C*E^2*u^4-480*B*C*E^2*u^3+1200*B*C*E^2*u^2+24*B*D^2*E*u^5+150*B*D^2*E*u^4+920*B*D^2*E*u^3+ 3300*B*D^2*E*u^2-96*B*D*E^2*u^5-80*B*D^3+1920*u*E^4+840*B*D*E^2*u^4+1680*B*D*E^2*u^2+9600*A^2*B*D*u+28560*A*B^2* E*u+3360*D*A^2*C*u+33600*B*A^2*E*u+960*A^2*B*C*u+7560*A*B^2*D*u+3360*A*B*D^2*u+6 720*A^2*E*D*u+9600*A*B*E^2*u+16800*A*C^2*E*u+16320*A*E^2*D*u+14400*A*E^2*C*u+420 0*A*C*D^2*u+360*A*C*B^2*u+12720*A*E*D^2*u+12960*C^2*E*D*u+12600*C*D^2*E*u+19320* B^2*C*E*u+13800*B^2*E*D*u-1470*A*B^2*D*u^4+8120*A*B^2*D*u^3-2100*A*B^2*D*u^2-840*A*B*D^2*u^5+4830*A*B*D^2*u^4-6600*A*B*D^2*u^3+12180*A*B*D^2*u^2-2688*A^2*E*D*u^5+13440*A^2*E*D*u^4-19200*A^2*E*D*u^3+16800*A^2*E*D*u^2-2688*A*B*E^2*u^5+13440*A*B*E^2*u^4-20160*A*B*E^2*u^3+16800*A*B*E^2*u^2+672*A*C^2*E*u^5-3360*A*C^2*E*u^4+11200*A*C^2*E*u^3-8640*A*C^2*E*u^2+192*A*E^2*D*u^5+1200*A*E^2*D*u^4-1600*A*E^2*D*u^3-5120*A^3*E-2800*B^2*A^2-2240*B*A^3-2240*A^3*C+25920*A*C*u*E*D+3780*D^3*u^2*A+2100*D^3*u^2*E+1800*D^3*u*A+4200*D^3*u *E+960*u^2*C^4-160*A*D*u^3*C^2+3360*u^2*C^3*D+3360*u^2*C^3*A-2992*A^3*D+6720*B*u*E^3)/(A+E+D+B+C)^(7/2)

Zugegeben, ein "wenig" länglich ist es schon. ;) Und die u-Terme (u, u^2, u^3, u^4, u^5) müssten noch ausgeklammert werden, damit die Darstellung schöner wird.

Mit einer vollständigen Lösung kann ich auch nicht dienen, Maple spuckt da auch was unschönes aus.

Neomi
2007-03-25, 21:01:51
... Taylorreihenentwicklung ...

Eigentlich hatte ich auf eine "handlichere" Möglichkeit gehofft. Irgendetwas, womit man eine Genauigkeit im Rahmen der Meßtoleranz bekommt, was aber nicht so lange dauert wie z.B. 64 Iterationen über die Tangentenlänge.

Nur so aus Interesse, ist "was unschönes" von Maple auch fehlerhaft, oder einfach nur extrem lang?

Spasstiger
2007-03-25, 21:36:22
Nur so aus Interesse, ist "was unschönes" von Maple auch fehlerhaft, oder einfach nur extrem lang?
Da kommt mehr oder weniger das raus, was auch dir von dir verlinkte Seite liefert.

Mit der Taylorreihenentwicklung hättest du halt letztendlich einfach nur Muliplikationen von u, u^2, u^3, u^4, etc. mit Faktoren, die wiederum Summen und Produkte deiner Paramter bilden. Mit fünfter Ordnung sollte eine Genauigkeit im Rahmen der Meßungenauigkeit schon ganz gut erreicht werden. Und vielleicht kann man durch geschicktes Zusammenfassen der Faktoren, die sich aus den Parametern zusammensetzen, auch noch ein paar Multiplikationen oder Additionen einsparen. Ich weiß ja nicht, wieviele dieser Operationen bei deinem Iterationsverfahren anfallen.

Wenn du die Taylorreihenentwicklung nur bis zur ersten Ordnung machst und dann integrierst, kommst du übrigens auf folgenden Ausdruck für das Integral:
http://www.250kb.de/u/070325/p/657f4f9a.png
Das wäre jetzt für die Entwicklung um den Punkt u=1, ist also nur dort einigermaßen brauchbar. Ich weiß ja nicht, in welchem Wertebereich dein u liegt. ;)

Wenn ich mit der Formel das Integral über u=0,9 bis u=1,1 ausrechne und für die Parameter A,B,C,D und E jeweils 1 einsetze, komme ich auf eine Abweichung von nur 0,17% vom richtigen Ergebniss. Für viele Anwendungen wäre das schon ausreichend genau.
Wenn natürlich deine Integralgrenzen deutlich weiter vom Entwicklungspunkt entfernt sind als nur meine +-10%, wird die Abweichung des Ergebnisses deutlich größer.

P.S. Ich lass Maple jetzt schon seit über 10 Minuten das Integral über die Grenzen von 0 bis 2 rechnen. Der Speichbedarf steigt stetig an und liegt schon bei mehr als 150 MB. Fieses Integral. Ein Core 2 Duo wäre jetzt auch ganz nützlich. Du kannst dir ja vorstellen, dass dieses Integral für Echtzeitrendering nicht mit Single-Precision-Genauigkeit oder gar Double-Precision-Genauigkeit gelöst werden sollte. ;)

/EDIT: Die Rechnung nähert sich den 20 Minuten, ein Ende ist nicht in Sicht.
/EDIT2: Schon 35 Minuten gerechnet und 275 MB Speicherbedarf. Ich brech die Rechnung jetzt ab. Hätte nicht gedacht, dass man mit so einem popeligen Integral den Rechner zur Weißglut treiben kann. ;)

Neomi
2007-03-25, 23:00:56
Holla, ganz schön rechenintensives Biest, das Integral. Das deutet doch sehr darauf hin, daß es keine exakte Lösung mit endlicher Länge gibt.

Der Wertebereich für die Integralgrenzen wird bei mir wohl überall zwischen 0 und 1 inklusive sein, also die normale Splineinterpolation, wobei ich eine Extrapolation mit Wertebereichen kleiner 0 und größer 1 nicht ausschließen kann. Single Precision ist definitiv ausreichend für das Anwendungsgebiet. Das war auch eher Neugier, ob es da noch einen besseren Weg gibt als iterative Näherungen.

Aber danke für die Mühe und Rechenzeit.

PS: ich hätte da noch einen Spline als Polynomial 5. Grades (Hermite mit Erweiterungen für C2-Stetigkeit zwischen einzelnen Splines), die Länge der Tangente läßt sich da auf diese Formel reduzieren:

sqrt(A*u^8 + B*u^7 + C*u^6 + D*u^5 + E*u^4 + F*u^3 + G*u^2 + H*u + I)

Wie stünden hier die Chancen für ein Integral? :D

Spasstiger
2007-03-25, 23:41:54
Wie stünden hier die Chancen für ein Integral? :D
Wird noch böser würde ich sagen. :biggrin:

Ich hätte da noch eine Idee, wie man den Rechenaufwand erträglich gestalten kann. Du stellst einfach mehrere Taylorentwicklungen für verschiedene Entwicklungspunkte zwischen 0 und 1 auf und bestimmst die zugehörigen Integrale. Wenn du die Taylorentwicklung bis zweiten Grades machst, hast du dann bei jeder Integralformel einen Bruch und eine Ordnung von u mehr gegenüber der Formel, die ich oben als Bildchen gepostet habe.
Dann hast du z.B. 10 Integralformeln mit jeweils vertretbarem Rechenaufwand und wählst jeweils nur die Integralformel aus, die zum u passt. Wenn du z.B. Integralformeln für die Entwicklungen um 1.0, 0.9. 0.8, 0.7, etc. aufgestellt hast und deine Grenzen liegen bei u=0.92 und u=0.37, dann wählst die einmal die Formel für 0.9 und einmal die Fomel für 0.4 aus. Auf diese Art solltest du die Rechenabweichungen im locker Zehntel-Promille-Bereich halten können (man liegt eher im hundertstel bis tausendstel Promille-Bereich, schätze ich). Ich weiß allerdings nicht, ob dir das als Genauigkeit genügt.

Wenn dir auch ein Promille Genauigkeit reicht, kannst du auch nur bis zur ersten Ordnung entwicklen und integrieren wie in der Formel oben im Bildchen. Das ganze musst du dann für zehn Entwicklungspunkte zwischen 0 und 1 machen (um eine Genauigkeit im Promillebereich zu erreichen).

Hier mal noch das Integral für die Entwicklung bis erster Ordnung um u=0.5:
http://www.250kb.de/u/070325/p/4d676500.png
Das sind zweimal je 11 Additionen/Subtraktionen, 12 Multiplikationen, eine Division und eine Quadratwurzel bis zum fertig gelösten Integral über einen bestimmten Bereich. Hinzu kommt noch zweimal die Routine zur Ermittlung der jeweils passenden Integralformel (Sprung) und eine Addition/Subtraktion (Ergebnis für obere Grenze minus Ergebnis für untere Grenze).
Das ist mit Sicherheit deutlich schneller als 64 Iterationen auszuführen.

Neomi
2007-03-26, 00:39:08
Ich weiß allerdings nicht, ob dir das als Genauigkeit genügt.

Wie schon angedeutet, wirklich nötig ist eine solche Genauigkeit für mich nicht. Das ist mehr eine Geschichte in Richtung "wenn schon, denn schon". Manchmal bastle ich mir Dinge, die ich wahrscheinlich nie brauchen werde, nur deshalb, weil es theoretisch möglich ist. ;)

Das sind zweimal je 11 Additionen/Subtraktionen, 12 Multiplikationen, eine Division und eine Quadratwurzel bis zum fertig gelösten Integral über einen bestimmten Bereich.

Einige Rechenschritte (inklusive der Wurzel, da kommt ja kein u drin vor) lassen sich auch rauskürzen bzw. einmal vorberechnen (und cachen) und mehrfach wiederverwenden. Ist also wohl eine ziemlich billige Methode.

Das ist mit Sicherheit deutlich schneller als 64 Iterationen auszuführen.

Jep, vor allem da jede Iteration noch mit einer Wurzel daherkommt. Und Wurzeln sind pöse. Die Wurzel ist ja eigentlich auch die einzige Sache, die mir das Integral versaut hat. Die Taylorreihen behalte ich auf jeden Fall im Hinterkopf.

Mit mehreren Reihen hat man leider ein Problem, wenn man andersrum rechnen will: man hat die Entfernung u=0 und will das dazu passende u haben. den passenden Bereich kennt man ja dann noch nicht. Kommt mir nur gerade in den Sinn.

Spasstiger
2007-03-26, 00:43:15
Die Taylorreihen behalte ich auf jeden Fall im Hinterkopf.
Taylorreihen sind generell eine nette Sache, wenn man nicht mit analytischen Methoden weiterkommt. :)

/EDIT: Hab noch bei Google gesucht und einen Thread in einem anderen Forum gefunden, der genau dasselbe behandelt:
http://www.uni-protokolle.de/foren/viewt/17931,0.html
Da wird von der Reihenentwicklung abgeraten, weil der Konvergenzbereich zu klein ist. Eine analytische Lösung wird aber auch dort nicht geliefert.