Die Variation der Elemente ist eine im 19. Jahrhundert entwickelte Methode zur genauen Bahnbestimmung von Himmelskörpern. Sie dient bis heute zur Modellierung so genannter Bahnstörungen.
In der Idealisierung des Zweikörperproblems verläuft die Umlaufbahn eines Planeten um die Sonne oder eines Mondes um einen Planeten exakt auf einer Keplerellipse. Voraussetzung hierfür ist, dass beide Körper kugelförmig sind, sich im Vakuum bewegen und keine weiteren Himmelskörper oder Kräfte wirken. Für die Berechnung solcher Keplerellipsen genügen sechs Bahnelemente und die drei Keplerschen Gesetze. Erstere 6 Zahlenwerte bleiben konstant - das heißt die Bahnellipse und ihre Ebene verändern sich nicht bezüglich des Zentralkörpers und des Fixsternhimmels.
De facto sind jedoch immer Bahnstörungen wirksam: dritte Körper, interplanetare Gase und Staub, Strahlungsdruck der Sonne, Atmosphären und Abplattung von Planeten usw. Diese störenden Kräfte verändern die 6 Bahnelemente langsam und verursachen zusätzliche Abweichungen von der Keplerbahn.
Lagrange und andere Astronomen entwickelten daher das Modell der oskulierenden Bahnen. Wenn die Umlaufbahn eines Himmelskörpers allzu variabel war, wird der momentan gültigen eine Ellipsenbahn angepasst, die sich allen Beobachtungen möglichst gut anschmiegt. So entsteht im Laufe der Zeit eine Folge oskulierender Bahnelemente, die stetig ineinander übergehen. Jeder dieser Elementensätze repräsentiert eine Bahn, auf welcher der Himmelskörper bei Aufhören der störenden Kräfte genau weiterfliegen würde.
Bei den Berechnungen werden näherungsweise mittels Differentialgleichungen die Elemente zu bestimmten Zeitpunkten ermittelt. Aus den errechneten Bahnen wird sodann die Hüllkurve bestimmt, welche dann die „theoretische“ Bahn des Objektes ergibt. Die zeitliche Veränderung der Bahnelemente kann säkular, periodisch und bis zu einem gewissen Grad auch unregelmäßig sein - je nach verursachender Kraft und Bahnelement. Für jeden der 6 Zahlenwerte lassen sich daher zeitabhängige Terme bestimmen, welche die Veränderung der Elemente charakterisieren.
Die Abweichungen der tatsächlichen Umlaufbahn von der knapp vorher gültigen oskulierenden Ellipse können als Funktion der störenden Krafte berechnet werden. Auf diese Art wurden erstmals 1846 Störungen der Uranusbahn modelliert, die zur Entdeckung des Planeten Neptun führten. Umgekehrt kann - zum Beispiel für Raumsonden - die Gravitation aller bekannten Körper berücksichtigt oder die Wirkung kurzer Korrekturmanöver berechnet werden.
Heutige Computer erlauben eine beliebig genaue Bahnbestimmung, wenn nur der Aufwand entsprechend hoch getrieben wird. Dabei bevorzugt man iterative Methoden und benützt zur sogenannten Bahnverbesserung die Ausgleichsrechnung bzw. die Kollokation. Nach einer ersten (genäherten) Bahnbestimmung aus wenigen Messungen (mindestens 3) berechnet man die Örter des Himmelskörpers zu den Zeitpunkten aller seiner Beobachtungen. Die Abweichungen von dieser Ephemeride sind die Eingangsgrößen der Bahnverbesserung. Sie bringt durch geeignete Variation der Elemente die berechneten Positionen in Einklang mit den Messungen. Dadurch und mittels Einbeziehung von Distanzmessungen mit Radar hat sich die Rechengenauigkeit im inneren Sonnensystem auf 1:10 Millionen und besser (Erdbahn auf km) erhöht.
Statt die Bahnelemente gezielt zu verändern, gibt es im Ergebnis gleichwertige Verfahren wie die Variation der geozentrischen Entfernung. Das Lösen der dabei entstehenden Differentialgleichungen erfolgt beispielsweise nach Runge-Kutta. Nach einer ähnlichen Methode arbeitet die numerische Integration des Jet Propulsion Laboratory (JPL-Programm DE200/LE200). Mit ihm werden alljährlich die Positionen aller großen Planeten und einiger Asteroiden für den Astronomical Almanac vorausberechnet.