Problem Summary

An ellipse with semiaxes \(a\) and \(b\) rolls on a straight line without slipping. Let \(C(a,b)\) denote the arc length of the trajectory traced by its center during one complete revolution. The problem asks for

$$C(1,4)+C(3,4).$$

The implementations do not search for an elementary antiderivative. Instead, they derive an exact speed formula for the center, reduce the geometry to a smooth one-dimensional integral, and evaluate that integral numerically to high precision.

Mathematical Approach

Work in the ellipse's own coordinate system and parameterize the boundary point that is instantaneously in contact with the ground by

$$r(t)=(x(t),y(t))=(a\cos t,b\sin t),\qquad 0\le t\le 2\pi.$$

As \(t\) runs once around the ellipse, the contact point visits the whole boundary once, so one full revolution of the rolling motion is captured by integrating over one full period of \(t\).

Step 1: Parameterize the ellipse and its local geometry

The first and second derivatives are

$$r'(t)=(-a\sin t,b\cos t),\qquad r''(t)=(-a\cos t,-b\sin t).$$

The distance from the center of the ellipse to the contact point is

$$\lVert r(t)\rVert=\sqrt{a^2\cos^2 t+b^2\sin^2 t}.$$

This distance varies with \(t\) unless \(a=b\), which is exactly why the center of a rolling ellipse does not move as simply as the center of a rolling circle.

Step 2: Turn rolling without slipping into a speed law

For pure rolling, the contact point has zero velocity relative to the ground at the instant of contact. Therefore that point is the instantaneous center of rotation of the rigid body.

If \(\theta(t)\) is the current orientation angle of the ellipse, then the center moves around the contact point with speed

$$v(t)=\left|\frac{d\theta}{dt}\right|\lVert r(t)\rVert.$$

So the geometric problem is reduced to finding the angular speed \(\theta'(t)\) associated with the changing tangent direction of the ellipse.

Step 3: Differentiate the tangent angle

For a smooth planar parametrized curve, the derivative of the tangent angle is

$$\frac{d\theta}{dt}=\frac{x'(t)y''(t)-y'(t)x''(t)}{x'(t)^2+y'(t)^2}.$$

Substituting the ellipse derivatives gives

$$x'(t)y''(t)-y'(t)x''(t)=ab,$$

$$x'(t)^2+y'(t)^2=a^2\sin^2 t+b^2\cos^2 t,$$

hence

$$\frac{d\theta}{dt}=\frac{ab}{a^2\sin^2 t+b^2\cos^2 t}.$$

The denominator is always positive for \(a,b>0\), so the angular speed stays positive and no sign ambiguity remains in the arc-length integral.

Step 4: Obtain the closed integrand for the center speed

Combining the previous two steps yields

$$v(t)=\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}.$$

Therefore the total center-path length over one revolution is

$$C(a,b)=\int_0^{2\pi}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

A useful sanity check is the circular case \(a=b=R\). Then the integrand becomes the constant \(R\), so

$$C(R,R)=\int_0^{2\pi}R\,dt=2\pi R,$$

which is exactly the expected path length for the center of a rolling circle.

Step 5: Use symmetry to integrate only one quarter-turn

The integrand depends only on \(\sin^2 t\) and \(\cos^2 t\). Consequently

$$v(t+\pi)=v(t),\qquad v(\pi-t)=v(t).$$

These symmetries imply

$$C(a,b)=4\int_0^{\pi/2}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

This is the exact form evaluated by the implementation, and it is numerically convenient because the interval is short and the integrand is smooth on \([0,\pi/2]\).

Step 6: Worked Example for \(C(2,4)\)

For \(a=2\) and \(b=4\), the quarter-turn integrand becomes

$$v(t)=\frac{4\sqrt{\cos^2 t+4\sin^2 t}}{\sin^2 t+4\cos^2 t}.$$

At the three standard Simpson nodes we get

$$v(0)=1,\qquad v\left(\frac{\pi}{4}\right)=\frac{4\sqrt{10}}{5},\qquad v\left(\frac{\pi}{2}\right)=8.$$

A single Simpson panel on \([0,\pi/2]\) gives the rough estimate

$$S=\frac{\pi/2}{6}\left(1+4\cdot\frac{4\sqrt{10}}{5}+8\right)=\frac{\pi}{12}\left(9+\frac{16\sqrt{10}}{5}\right)\approx 5.0056.$$

Multiplying by \(4\) would give about \(20.0224\), which is not yet accurate enough. After adaptive refinement the numerical integral converges to

$$C(2,4)\approx 21.38816906,$$

showing why an adaptive method is used instead of a single coarse panel.

How the Code Works

The C++, Python, and Java implementations all follow the same algorithm. First they evaluate the closed-form speed formula above. Next they apply Simpson's rule on the quarter interval \([0,\pi/2]\) using the endpoint values and the midpoint value.

They then recurse adaptively. For a current interval, the implementation evaluates the speed at the two quarter points, compares the sum of the left and right Simpson estimates with the Simpson estimate on the whole interval, and interprets the difference as an error signal.

If the estimated local error already satisfies

$$|\Delta|\le 15\varepsilon,$$

or if the recursion limit has been reached, the code accepts the interval and applies the standard Richardson correction

$$S_{\text{refined}}=S_{\text{left}}+S_{\text{right}}+\frac{\Delta}{15}.$$

Otherwise it splits the interval in half and repeats the same process on both subintervals. The requested integral is obtained by multiplying the converged quarter-interval value by \(4\). Finally the implementation evaluates \(C(1,4)\) and \(C(3,4)\), adds them, and prints the result to eight decimal places. The numerical parameters are \(\varepsilon=10^{-13}\) and a maximum recursion depth of \(30\).

Complexity Analysis

Adaptive Simpson integration has no fixed iteration count in advance, so it is best analyzed in terms of the number of accepted subintervals. If the recursion finishes with \(N\) accepted leaves, the total time is \(O(N)\) function evaluations up to constant factors, because each subdivision adds only a constant number of new samples. The memory usage is \(O(D)\), where \(D\) is the recursion depth; in the implementation \(D\le 30\).

For this problem the integrand is smooth on \([0,\pi/2]\), so \(N\) stays modest in practice. The adaptive strategy is efficient because it refines only where the shape of the integrand demands it, instead of spending the same resolution everywhere.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=525
  2. Ellipse: Wikipedia — Ellipse
  3. Arc length: Wikipedia — Arc length
  4. Adaptive Simpson's method: Wikipedia — Adaptive Simpson's method

Problemzusammenfassung

Eine Ellipse mit Halbachsen \(a\) und \(b\) rollt ohne Gleiten auf einer Geraden. Mit \(C(a,b)\) bezeichnen wir die Bogenlänge der Bahn, die ihr Mittelpunkt während einer vollständigen Umdrehung zurücklegt. Gesucht ist

$$C(1,4)+C(3,4).$$

Die Implementierungen suchen nicht nach einer elementaren Stammfunktion. Stattdessen wird zuerst eine exakte Geschwindigkeitsformel für den Mittelpunkt hergeleitet, die Geometrie auf ein glattes eindimensionales Integral reduziert und dieses Integral anschließend numerisch mit hoher Genauigkeit ausgewertet.

Mathematischer Ansatz

Wir arbeiten im körperfesten Koordinatensystem der Ellipse und parametrisieren den Randpunkt, der gerade den Boden berührt, durch

$$r(t)=(x(t),y(t))=(a\cos t,b\sin t),\qquad 0\le t\le 2\pi.$$

Wenn \(t\) die Ellipse einmal umläuft, wird der gesamte Rand genau einmal durchlaufen. Daher wird eine vollständige Rollbewegung durch die Integration über eine volle Periode von \(t\) beschrieben.

Schritt 1: Parametrisiere die Ellipse und ihre lokale Geometrie

Die erste und zweite Ableitung lauten

$$r'(t)=(-a\sin t,b\cos t),\qquad r''(t)=(-a\cos t,-b\sin t).$$

Der Abstand vom Mittelpunkt der Ellipse zum Berührpunkt ist

$$\lVert r(t)\rVert=\sqrt{a^2\cos^2 t+b^2\sin^2 t}.$$

Dieser Abstand hängt von \(t\) ab, solange \(a\ne b\). Genau deshalb verhält sich der Mittelpunkt einer rollenden Ellipse wesentlich komplizierter als der Mittelpunkt eines rollenden Kreises.

Schritt 2: Wandle „ohne Gleiten“ in ein Geschwindigkeitsgesetz um

Bei reinem Rollen hat der Berührpunkt im Kontaktmoment relativ zum Boden Geschwindigkeit \(0\). Also ist dieser Punkt das momentane Drehzentrum des starren Körpers.

Bezeichnet \(\theta(t)\) den aktuellen Orientierungswinkel der Ellipse, dann bewegt sich der Mittelpunkt um den Berührpunkt mit der Geschwindigkeit

$$v(t)=\left|\frac{d\theta}{dt}\right|\lVert r(t)\rVert.$$

Damit ist das geometrische Problem darauf reduziert, die Winkelgeschwindigkeit \(\theta'(t)\) zu bestimmen, die aus der Änderung der Tangentenrichtung entsteht.

Schritt 3: Leite den Tangentenwinkel ab

Für eine glatte ebene Parametrisierung gilt für die Ableitung des Tangentenwinkels

$$\frac{d\theta}{dt}=\frac{x'(t)y''(t)-y'(t)x''(t)}{x'(t)^2+y'(t)^2}.$$

Für die Ellipse erhält man damit

$$x'(t)y''(t)-y'(t)x''(t)=ab,$$

$$x'(t)^2+y'(t)^2=a^2\sin^2 t+b^2\cos^2 t,$$

also

$$\frac{d\theta}{dt}=\frac{ab}{a^2\sin^2 t+b^2\cos^2 t}.$$

Der Nenner ist für \(a,b>0\) stets positiv. Daher bleibt auch die Winkelgeschwindigkeit positiv, und im Bogenlängenintegral gibt es kein Vorzeichenproblem.

Schritt 4: Erhalte den geschlossenen Integranden für die Mittelpunktgeschwindigkeit

Aus den beiden vorherigen Schritten folgt

$$v(t)=\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}.$$

Damit ist die gesamte Bahnlänge des Mittelpunkts nach einer vollen Umdrehung

$$C(a,b)=\int_0^{2\pi}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

Ein nützlicher Plausibilitätscheck ist der Kreisfall \(a=b=R\). Dann ist der Integrand konstant gleich \(R\), also

$$C(R,R)=\int_0^{2\pi}R\,dt=2\pi R,$$

genau die erwartete Weglänge des Mittelpunkts eines rollenden Kreises.

Schritt 5: Nutze Symmetrie und integriere nur über ein Viertel

Der Integrand hängt nur von \(\sin^2 t\) und \(\cos^2 t\) ab. Deshalb gilt

$$v(t+\pi)=v(t),\qquad v(\pi-t)=v(t).$$

Daraus folgt

$$C(a,b)=4\int_0^{\pi/2}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

Genau diese Form wird in der Implementierung ausgewertet. Numerisch ist sie günstig, weil das Integrationsintervall kurz ist und der Integrand auf \([0,\pi/2]\) glatt bleibt.

Schritt 6: Durchgerechnetes Beispiel für \(C(2,4)\)

Für \(a=2\) und \(b=4\) wird der Integrand des Viertelintervalls zu

$$v(t)=\frac{4\sqrt{\cos^2 t+4\sin^2 t}}{\sin^2 t+4\cos^2 t}.$$

An den drei klassischen Simpson-Stützstellen erhält man

$$v(0)=1,\qquad v\left(\frac{\pi}{4}\right)=\frac{4\sqrt{10}}{5},\qquad v\left(\frac{\pi}{2}\right)=8.$$

Ein einziges Simpson-Panel auf \([0,\pi/2]\) liefert die grobe Näherung

$$S=\frac{\pi/2}{6}\left(1+4\cdot\frac{4\sqrt{10}}{5}+8\right)=\frac{\pi}{12}\left(9+\frac{16\sqrt{10}}{5}\right)\approx 5.0056.$$

Mit dem Faktor \(4\) ergäbe das nur etwa \(20.0224\) und ist damit noch nicht präzise genug. Nach adaptiver Verfeinerung konvergiert das numerische Integral dagegen zu

$$C(2,4)\approx 21.38816906,$$

und genau daran sieht man, weshalb ein adaptives Verfahren anstelle eines einzelnen groben Panels verwendet wird.

Wie der Code arbeitet

Die Implementierungen in C++, Python und Java verwenden denselben Ablauf. Zuerst werten sie die obige geschlossene Geschwindigkeitsformel aus. Danach wird auf dem Viertelintervall \([0,\pi/2]\) die Simpson-Regel mit den beiden Randwerten und dem Mittelpunkt angewandt.

Anschließend folgt eine adaptive Rekursion. Für ein aktuelles Intervall werden die Geschwindigkeiten an den beiden Viertelpunkten ausgewertet, die Summe der linken und rechten Simpson-Näherung mit der Simpson-Näherung des ganzen Intervalls verglichen und die Differenz als Fehlersignal interpretiert.

Falls der geschätzte lokale Fehler bereits

$$|\Delta|\le 15\varepsilon$$

erfüllt oder die Rekursionstiefe ausgeschöpft ist, akzeptiert der Code das Intervall und verwendet die übliche Richardson-Korrektur

$$S_{\text{refined}}=S_{\text{left}}+S_{\text{right}}+\frac{\Delta}{15}.$$

Andernfalls wird das Intervall halbiert und derselbe Prozess auf beide Teilintervalle angewandt. Das gesuchte Integral entsteht aus dem konvergierten Viertelwert durch Multiplikation mit \(4\). Danach berechnet die Implementierung \(C(1,4)\) und \(C(3,4)\), addiert beide Werte und gibt das Ergebnis mit acht Nachkommastellen aus. Verwendet werden \(\varepsilon=10^{-13}\) und eine maximale Rekursionstiefe von \(30\).

Komplexitätsanalyse

Für adaptive Simpson-Integration gibt es keine im Voraus feste Iterationszahl. Am sinnvollsten analysiert man daher die Methode über die Anzahl akzeptierter Teilintervalle. Wenn die Rekursion mit \(N\) akzeptierten Blättern endet, beträgt die Laufzeit bis auf konstante Faktoren \(O(N)\) Funktionsauswertungen, denn jede Verfeinerung fügt nur konstant viele neue Stützstellen hinzu. Der Speicherverbrauch ist \(O(D)\), wobei \(D\) die Rekursionstiefe bezeichnet; in der Implementierung gilt \(D\le 30\).

In dieser Aufgabe ist der Integrand auf \([0,\pi/2]\) glatt, daher bleibt \(N\) in der Praxis klein. Die adaptive Strategie ist effizient, weil sie nur dort weiter verfeinert, wo die Form des Integranden es verlangt, statt überall dieselbe Auflösung zu erzwingen.

Fußnoten und Quellen

  1. Problemseite: https://projecteuler.net/problem=525
  2. Ellipse: Wikipedia — Ellipse
  3. Bogenlänge: Wikipedia — Arc length
  4. Adaptive Simpson-Methode: Wikipedia — Adaptive Simpson's method

Problem Özeti

Yarı eksenleri \(a\) ve \(b\) olan bir elips, bir doğru üzerinde kaymadan yuvarlanıyor. \(C(a,b)\), merkezin bir tam dönüş boyunca izlediği yolun yay uzunluğu olsun. İstenen değer

$$C(1,4)+C(3,4).$$

Uygulamalar kapalı biçimli elementer bir asal integral aramıyor. Bunun yerine önce merkez hızının tam formülü çıkarılıyor, geometri düzgün bir tek boyutlu integrale indirgeniyor ve bu integral yüksek doğrulukla sayısal olarak hesaplanıyor.

Matematiksel Yaklaşım

Elipsin kendi koordinat sisteminde çalışalım ve o anda yere temas eden sınır noktasını

$$r(t)=(x(t),y(t))=(a\cos t,b\sin t),\qquad 0\le t\le 2\pi$$

ile parametrize edelim. \(t\) bir kez tüm elipsi dolaştığında temas noktası da sınırın tamamını bir kez gezer; dolayısıyla tam bir yuvarlanma hareketi \(t\)'nin bir tam periyodu üzerinde entegrasyon yapılarak yakalanır.

Adım 1: Elipsi ve yerel geometrisini parametrize et

Birinci ve ikinci türevler

$$r'(t)=(-a\sin t,b\cos t),\qquad r''(t)=(-a\cos t,-b\sin t)$$

şeklindedir. Elips merkezinden temas noktasına olan uzaklık ise

$$\lVert r(t)\rVert=\sqrt{a^2\cos^2 t+b^2\sin^2 t}$$

olur. \(a=b\) olmadıkça bu uzaklık sabit değildir; yuvarlanan elipsin merkez hareketinin daire durumundan daha karmaşık olmasının nedeni budur.

Adım 2: Kaymasız yuvarlanmayı hız yasasına çevir

Saf yuvarlanmada temas noktası, temas anında yere göre sıfır hıza sahiptir. Bu yüzden o nokta rijit cismin anlık dönme merkezidir.

\(\theta(t)\) elipsin o andaki yönelim açısı olsun. O zaman merkezin temas noktasına göre hızı

$$v(t)=\left|\frac{d\theta}{dt}\right|\lVert r(t)\rVert$$

olur. Böylece geometrik problem, teğet yönündeki değişimden doğan açısal hız \(\theta'(t)\)'yi bulma problemine indirgenir.

Adım 3: Teğet açısının türevini hesapla

Düzlemde düzgün parametrize edilmiş bir eğri için teğet açısının türevi

$$\frac{d\theta}{dt}=\frac{x'(t)y''(t)-y'(t)x''(t)}{x'(t)^2+y'(t)^2}$$

eşitliğiyle verilir. Elips için yerine koyarsak

$$x'(t)y''(t)-y'(t)x''(t)=ab,$$

$$x'(t)^2+y'(t)^2=a^2\sin^2 t+b^2\cos^2 t$$

ve dolayısıyla

$$\frac{d\theta}{dt}=\frac{ab}{a^2\sin^2 t+b^2\cos^2 t}$$

elde edilir. \(a,b>0\) için payda her zaman pozitiftir; bu yüzden açısal hız da pozitiftir ve yay uzunluğu integralinde işaret belirsizliği kalmaz.

Adım 4: Merkez hızı için kapalı biçimli integrandı elde et

Önceki iki adımı birleştirince

$$v(t)=\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}$$

bulunur. Böylece bir tam dönüşte merkez yolunun uzunluğu

$$C(a,b)=\int_0^{2\pi}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt$$

olur. Faydalı bir sağlamlık kontrolü daire durumudur: \(a=b=R\) ise integrand sabit olarak \(R\)'ye iner ve

$$C(R,R)=\int_0^{2\pi}R\,dt=2\pi R$$

elde edilir; bu da yuvarlanan bir dairenin merkezi için beklenen yol uzunluğudur.

Adım 5: Simetriyi kullan ve sadece çeyrek turu integre et

İntegrand yalnızca \(\sin^2 t\) ve \(\cos^2 t\) içerir. Bu nedenle

$$v(t+\pi)=v(t),\qquad v(\pi-t)=v(t)$$

eşitlikleri geçerlidir. Buradan

$$C(a,b)=4\int_0^{\pi/2}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt$$

sonucuna ulaşılır. Uygulamanın değerlendirdiği tam formül budur; aralık kısa ve integrand \([0,\pi/2]\) üzerinde düzgün olduğu için sayısal olarak da elverişlidir.

Adım 6: \(C(2,4)\) için çözümlü örnek

\(a=2\) ve \(b=4\) için çeyrek tur integrandı

$$v(t)=\frac{4\sqrt{\cos^2 t+4\sin^2 t}}{\sin^2 t+4\cos^2 t}$$

olur. Üç klasik Simpson düğümünde

$$v(0)=1,\qquad v\left(\frac{\pi}{4}\right)=\frac{4\sqrt{10}}{5},\qquad v\left(\frac{\pi}{2}\right)=8$$

elde edilir. \([0,\pi/2]\) üzerinde tek Simpson paneli kaba tahmin olarak

$$S=\frac{\pi/2}{6}\left(1+4\cdot\frac{4\sqrt{10}}{5}+8\right)=\frac{\pi}{12}\left(9+\frac{16\sqrt{10}}{5}\right)\approx 5.0056$$

verir. Bunu \(4\) ile çarpmak yaklaşık \(20.0224\) eder ve henüz yeterince doğru değildir. Adaptif iyileştirme sonrasında sayısal integral

$$C(2,4)\approx 21.38816906$$

değerine yakınsar; tek bir kaba panel yerine neden adaptif yöntemin seçildiği böylece açıkça görülür.

Kod Nasıl Çalışır

C++, Python ve Java uygulamaları aynı algoritmayı izler. Önce yukarıdaki kapalı biçimli hız formülü değerlendirilir. Ardından \([0,\pi/2]\) çeyrek aralığında uç noktalar ve orta nokta kullanılarak Simpson kuralı uygulanır.

Daha sonra adaptif özyineleme başlar. Mevcut bir aralık için uygulama iki çeyrek noktadaki hızları hesaplar, sol ve sağ Simpson tahminlerinin toplamını tüm aralığın Simpson tahminiyle karşılaştırır ve farkı hata sinyali olarak kullanır.

Eğer yerel hata tahmini

$$|\Delta|\le 15\varepsilon$$

koşulunu sağlıyorsa ya da özyineleme sınırına gelinmişse, aralık kabul edilir ve standart Richardson düzeltmesi

$$S_{\text{refined}}=S_{\text{left}}+S_{\text{right}}+\frac{\Delta}{15}$$

uygulanır. Aksi halde aralık ikiye bölünür ve aynı süreç iki alt aralıkta tekrarlanır. İstenen integral, yakınsamış çeyrek-aralık sonucunun \(4\) ile çarpılmasıyla elde edilir. Son olarak uygulama \(C(1,4)\) ve \(C(3,4)\) değerlerini hesaplar, toplar ve sonucu sekiz ondalık basamakla yazar. Sayısal parametreler \(\varepsilon=10^{-13}\) ve en fazla \(30\) özyineleme derinliğidir.

Karmaşıklık Analizi

Adaptif Simpson entegrasyonunda yineleme sayısı önceden sabit değildir; bu yüzden en doğal analiz, kabul edilen alt aralık sayısı üzerinden yapılır. Özyineleme \(N\) adet kabul edilmiş yaprak aralıkla biterse toplam süre sabit çarpanlar dışında \(O(N)\) fonksiyon değerlendirmesidir; çünkü her bölme yalnızca sabit sayıda yeni örnek ekler. Bellek kullanımı \(O(D)\)'dir; burada \(D\) özyineleme derinliğidir ve uygulamada \(D\le 30\).

Bu problemde integrand \([0,\pi/2]\) üzerinde düzgündür, dolayısıyla pratikte \(N\) küçüktür. Adaptif strateji yalnızca integrandın biçiminin gerçekten gerektirdiği bölgelerde inceltme yaptığı için etkilidir.

Dipnotlar ve Kaynakça

  1. Problem sayfası: https://projecteuler.net/problem=525
  2. Ellipse: Wikipedia — Ellipse
  3. Yay uzunluğu: Wikipedia — Arc length
  4. Adaptif Simpson yöntemi: Wikipedia — Adaptive Simpson's method

Resumen del Problema

Una elipse de semiejes \(a\) y \(b\) rueda sobre una recta sin deslizamiento. Sea \(C(a,b)\) la longitud de arco de la trayectoria recorrida por su centro durante una revolución completa. Hay que calcular

$$C(1,4)+C(3,4).$$

Las implementaciones no buscan una primitiva elemental cerrada. En cambio, primero derivan una fórmula exacta para la velocidad del centro, reducen la geometría a una integral unidimensional suave y luego evalúan esa integral numéricamente con alta precisión.

Enfoque Matemático

Trabajamos en el sistema de coordenadas propio de la elipse y parametrizamos el punto del borde que en ese instante está en contacto con el suelo mediante

$$r(t)=(x(t),y(t))=(a\cos t,b\sin t),\qquad 0\le t\le 2\pi.$$

Cuando \(t\) recorre una vez la elipse, el punto de contacto visita toda la frontera una sola vez, así que una revolución completa del movimiento de rodadura queda descrita integrando sobre un período completo de \(t\).

Paso 1: Parametrizar la elipse y su geometría local

La primera y la segunda derivada son

$$r'(t)=(-a\sin t,b\cos t),\qquad r''(t)=(-a\cos t,-b\sin t).$$

La distancia desde el centro de la elipse hasta el punto de contacto es

$$\lVert r(t)\rVert=\sqrt{a^2\cos^2 t+b^2\sin^2 t}.$$

Esa distancia no es constante salvo cuando \(a=b\). Precisamente por eso el centro de una elipse rodante no se mueve de forma tan simple como el centro de un círculo rodante.

Paso 2: Convertir la rodadura sin deslizamiento en una ley de velocidad

En rodadura pura, el punto de contacto tiene velocidad cero respecto al suelo en el instante del contacto. Por tanto, ese punto es el centro instantáneo de rotación del sólido rígido.

Si \(\theta(t)\) es el ángulo de orientación actual de la elipse, entonces el centro se mueve alrededor del punto de contacto con velocidad

$$v(t)=\left|\frac{d\theta}{dt}\right|\lVert r(t)\rVert.$$

Así, el problema geométrico queda reducido a hallar la velocidad angular \(\theta'(t)\) producida por el cambio de dirección de la tangente de la elipse.

Paso 3: Derivar el ángulo de la tangente

Para una curva plana suave parametrizada, la derivada del ángulo de la tangente viene dada por

$$\frac{d\theta}{dt}=\frac{x'(t)y''(t)-y'(t)x''(t)}{x'(t)^2+y'(t)^2}.$$

Sustituyendo las derivadas de la elipse se obtiene

$$x'(t)y''(t)-y'(t)x''(t)=ab,$$

$$x'(t)^2+y'(t)^2=a^2\sin^2 t+b^2\cos^2 t,$$

y por lo tanto

$$\frac{d\theta}{dt}=\frac{ab}{a^2\sin^2 t+b^2\cos^2 t}.$$

El denominador es siempre positivo para \(a,b>0\), así que la velocidad angular también es positiva y no aparece ambigüedad de signo en la integral de longitud de arco.

Paso 4: Obtener el integrando cerrado para la velocidad del centro

Combinando los dos pasos anteriores resulta

$$v(t)=\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}.$$

En consecuencia, la longitud total de la trayectoria del centro en una revolución es

$$C(a,b)=\int_0^{2\pi}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

Una comprobación útil es el caso circular \(a=b=R\). Entonces el integrando se reduce a la constante \(R\), de modo que

$$C(R,R)=\int_0^{2\pi}R\,dt=2\pi R,$$

que es exactamente la distancia esperada para el centro de un círculo que rueda.

Paso 5: Usar simetría e integrar solo un cuarto de vuelta

El integrando depende solo de \(\sin^2 t\) y \(\cos^2 t\). Por eso se cumple

$$v(t+\pi)=v(t),\qquad v(\pi-t)=v(t).$$

Estas simetrías implican

$$C(a,b)=4\int_0^{\pi/2}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

Ésta es exactamente la forma que evalúa la implementación, y resulta cómoda numéricamente porque el intervalo es corto y el integrando es suave en \([0,\pi/2]\).

Paso 6: Ejemplo trabajado para \(C(2,4)\)

Para \(a=2\) y \(b=4\), el integrando en el cuarto de vuelta se convierte en

$$v(t)=\frac{4\sqrt{\cos^2 t+4\sin^2 t}}{\sin^2 t+4\cos^2 t}.$$

En los tres nodos clásicos de Simpson obtenemos

$$v(0)=1,\qquad v\left(\frac{\pi}{4}\right)=\frac{4\sqrt{10}}{5},\qquad v\left(\frac{\pi}{2}\right)=8.$$

Un único panel de Simpson en \([0,\pi/2]\) produce la estimación gruesa

$$S=\frac{\pi/2}{6}\left(1+4\cdot\frac{4\sqrt{10}}{5}+8\right)=\frac{\pi}{12}\left(9+\frac{16\sqrt{10}}{5}\right)\approx 5.0056.$$

Multiplicar por \(4\) da aproximadamente \(20.0224\), que todavía no es suficientemente preciso. Tras el refinamiento adaptativo, la integral numérica converge a

$$C(2,4)\approx 21.38816906,$$

lo que muestra por qué se usa un método adaptativo en lugar de un panel grueso único.

Cómo Funciona el Código

Las implementaciones en C++, Python y Java siguen exactamente el mismo esquema. Primero evalúan la fórmula cerrada de la velocidad. Después aplican la regla de Simpson en el intervalo de cuarto de vuelta \([0,\pi/2]\) usando los dos extremos y el punto medio.

Luego entra en juego la recursión adaptativa. Para un intervalo dado, la implementación evalúa la velocidad en los dos cuartos interiores, compara la suma de las aproximaciones de Simpson izquierda y derecha con la aproximación de Simpson del intervalo completo e interpreta la diferencia como una señal de error.

Si el error local estimado ya satisface

$$|\Delta|\le 15\varepsilon,$$

o si se alcanzó el límite de recursión, el código acepta el intervalo y aplica la corrección estándar de Richardson

$$S_{\text{refined}}=S_{\text{left}}+S_{\text{right}}+\frac{\Delta}{15}.$$

En caso contrario, divide el intervalo en dos y repite el proceso en ambos subintervalos. La integral buscada se obtiene multiplicando por \(4\) el valor convergido en el cuarto de intervalo. Por último, la implementación calcula \(C(1,4)\) y \(C(3,4)\), suma ambos valores y muestra el resultado con ocho cifras decimales. Los parámetros numéricos elegidos son \(\varepsilon=10^{-13}\) y profundidad máxima de recursión \(30\).

Análisis de Complejidad

La integración adaptativa de Simpson no tiene un número fijo de iteraciones de antemano, así que conviene analizarla en función del número de subintervalos aceptados. Si la recursión termina con \(N\) hojas aceptadas, el tiempo total es \(O(N)\) evaluaciones de la función salvo factores constantes, porque cada subdivisión añade solo un número constante de nuevas muestras. El uso de memoria es \(O(D)\), donde \(D\) es la profundidad de la recursión; en la implementación se cumple \(D\le 30\).

En este problema el integrando es suave en \([0,\pi/2]\), así que \(N\) se mantiene pequeño en la práctica. La estrategia adaptativa es eficiente porque refina solo donde la forma del integrando realmente lo exige, en vez de imponer la misma resolución en todo el intervalo.

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=525
  2. Elipse: Wikipedia — Ellipse
  3. Longitud de arco: Wikipedia — Arc length
  4. Método de Simpson adaptativo: Wikipedia — Adaptive Simpson's method

问题概述

半轴分别为 \(a\) 和 \(b\) 的椭圆在一条直线上做无滑动滚动。记 \(C(a,b)\) 为它的中心在一次完整转动过程中所走轨迹的弧长。题目要求计算

$$C(1,4)+C(3,4).$$

实现并没有去寻找一个初等闭式原函数,而是先把几何问题化成中心速度的精确表达式,再把整个问题归结为一个光滑的一维定积分,最后用高精度自适应数值积分求值。

数学方法

在椭圆自身的坐标系中,设当前与地面接触的边界点参数化为

$$r(t)=(x(t),y(t))=(a\cos t,b\sin t),\qquad 0\le t\le 2\pi.$$

当 \(t\) 从 \(0\) 走到 \(2\pi\) 时,这个接触点会把整个椭圆边界恰好走一遍,因此一次完整的滚动也就可以通过对 \(t\) 的一个完整周期积分来描述。

步骤 1:参数化椭圆并写出局部几何

一阶和二阶导数为

$$r'(t)=(-a\sin t,b\cos t),\qquad r''(t)=(-a\cos t,-b\sin t).$$

椭圆中心到接触点的距离是

$$\lVert r(t)\rVert=\sqrt{a^2\cos^2 t+b^2\sin^2 t}.$$

除非 \(a=b\),否则这个距离会随 \(t\) 改变。也正因为如此,滚动椭圆的中心轨迹不会像滚动圆的中心那样简单。

步骤 2:把无滑动滚动转成速度公式

纯滚动意味着接触点在接触瞬间相对于地面的速度为零,所以该点就是刚体的瞬时转动中心。

若 \(\theta(t)\) 表示椭圆此时的朝向角,那么中心绕接触点运动的速度满足

$$v(t)=\left|\frac{d\theta}{dt}\right|\lVert r(t)\rVert.$$

于是问题就转化成:求出由于椭圆切线方向变化所产生的角速度 \(\theta'(t)\)。

步骤 3:求切线方向角的导数

对于一条光滑的平面参数曲线,切线方向角的导数公式为

$$\frac{d\theta}{dt}=\frac{x'(t)y''(t)-y'(t)x''(t)}{x'(t)^2+y'(t)^2}.$$

把椭圆的导数代入后得到

$$x'(t)y''(t)-y'(t)x''(t)=ab,$$

$$x'(t)^2+y'(t)^2=a^2\sin^2 t+b^2\cos^2 t,$$

因此

$$\frac{d\theta}{dt}=\frac{ab}{a^2\sin^2 t+b^2\cos^2 t}.$$

当 \(a,b>0\) 时分母始终为正,所以角速度始终为正,弧长积分中也就不再有符号上的歧义。

步骤 4:得到中心速度的闭式被积函数

把前两步合并起来,中心速度就变成

$$v(t)=\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}.$$

于是一次完整转动中中心路径的长度为

$$C(a,b)=\int_0^{2\pi}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

这里有一个很好的自检:如果 \(a=b=R\),也就是退化成圆,那么被积函数就恒等于 \(R\),从而

$$C(R,R)=\int_0^{2\pi}R\,dt=2\pi R,$$

这正是滚动圆的中心应当走出的距离。

步骤 5:利用对称性只积分四分之一圈

由于被积函数只含有 \(\sin^2 t\) 和 \(\cos^2 t\),所以它满足

$$v(t+\pi)=v(t),\qquad v(\pi-t)=v(t).$$

于是可把整圈积分化简为

$$C(a,b)=4\int_0^{\pi/2}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

实现里实际计算的正是这个式子。这样做的好处是积分区间短,而且 \([0,\pi/2]\) 上的被积函数非常平滑,适合做高精度数值积分。

步骤 6:\(C(2,4)\) 的算例

取 \(a=2\)、\(b=4\) 时,四分之一圈上的被积函数化为

$$v(t)=\frac{4\sqrt{\cos^2 t+4\sin^2 t}}{\sin^2 t+4\cos^2 t}.$$

在 Simpson 公式最常用的三个节点上,有

$$v(0)=1,\qquad v\left(\frac{\pi}{4}\right)=\frac{4\sqrt{10}}{5},\qquad v\left(\frac{\pi}{2}\right)=8.$$

如果只在 \([0,\pi/2]\) 上做一个 Simpson 面板,那么得到的粗略估计是

$$S=\frac{\pi/2}{6}\left(1+4\cdot\frac{4\sqrt{10}}{5}+8\right)=\frac{\pi}{12}\left(9+\frac{16\sqrt{10}}{5}\right)\approx 5.0056.$$

再乘上外面的 \(4\) 后只有大约 \(20.0224\),精度显然不够。继续做自适应细分以后,数值积分会收敛到

$$C(2,4)\approx 21.38816906,$$

这也说明了为什么实现要用自适应方法,而不是只用一个很粗的 Simpson 面板。

代码如何工作

C++、Python 和 Java 三个实现采用的是同一套流程。第一步是计算上面的闭式速度公式。第二步是在四分之一圈区间 \([0,\pi/2]\) 上用端点和中点做一次 Simpson 近似。

接着进入自适应递归。对当前区间,程序会额外计算两个四分点处的速度值,把左右两个半区间的 Simpson 结果相加,再与整个区间的 Simpson 结果比较,把差值当作局部误差信号。

如果估计误差已经满足

$$|\Delta|\le 15\varepsilon,$$

或者递归深度已经到达上限,那么该区间就被接受,并使用标准的 Richardson 修正

$$S_{\text{refined}}=S_{\text{left}}+S_{\text{right}}+\frac{\Delta}{15}.$$

否则就把区间二分,并在两个子区间上重复同样的过程。收敛后的四分之一圈积分值再乘以 \(4\),就得到 \(C(a,b)\)。最后,程序分别计算 \(C(1,4)\) 和 \(C(3,4)\),把它们相加,并以八位小数输出。实现中使用的数值参数是 \(\varepsilon=10^{-13}\),递归深度上限为 \(30\)。

复杂度分析

自适应 Simpson 积分没有预先固定的迭代次数,因此最自然的分析方式是看最终接受了多少个子区间。若递归结束时共有 \(N\) 个被接受的叶子区间,那么总时间复杂度在常数因子意义下是 \(O(N)\) 次函数求值,因为每次细分只会增加常数个新的采样点。空间复杂度为 \(O(D)\),其中 \(D\) 是递归深度;在实现里有 \(D\le 30\)。

本题中的被积函数在 \([0,\pi/2]\) 上很光滑,所以实际需要的 \(N\) 并不大。自适应策略之所以高效,是因为它只在真正需要更高分辨率的局部区域继续细分,而不是在整个区间上平均地浪费采样点。

注释与参考资料

  1. 题目页面: https://projecteuler.net/problem=525
  2. Ellipse: Wikipedia — Ellipse
  3. 弧长: Wikipedia — Arc length
  4. 自适应 Simpson 方法: Wikipedia — Adaptive Simpson's method

Краткое описание задачи

Эллипс с полуосями \(a\) и \(b\) катится по прямой без скольжения. Обозначим через \(C(a,b)\) длину траектории, которую проходит его центр за один полный оборот. Нужно вычислить

$$C(1,4)+C(3,4).$$

Реализации не ищут элементарную первообразную в замкнутом виде. Вместо этого сначала выводится точная формула для скорости центра, затем геометрия сводится к гладкому одномерному интегралу, и этот интеграл вычисляется численно с высокой точностью.

Математический подход

Будем работать в собственной системе координат эллипса и параметризуем точку границы, которая в данный момент касается земли, так:

$$r(t)=(x(t),y(t))=(a\cos t,b\sin t),\qquad 0\le t\le 2\pi.$$

Когда \(t\) один раз проходит всю окружность параметра, точка контакта обходит всю границу эллипса ровно один раз. Поэтому один полный цикл качения описывается интегрированием по полному периоду параметра \(t\).

Шаг 1: Параметризовать эллипс и выписать локальную геометрию

Первая и вторая производные имеют вид

$$r'(t)=(-a\sin t,b\cos t),\qquad r''(t)=(-a\cos t,-b\sin t).$$

Расстояние от центра эллипса до точки контакта равно

$$\lVert r(t)\rVert=\sqrt{a^2\cos^2 t+b^2\sin^2 t}.$$

Если только \(a\ne b\), это расстояние меняется вместе с \(t\). Именно поэтому движение центра катящегося эллипса существенно сложнее, чем движение центра катящегося круга.

Шаг 2: Преобразовать условие качения без скольжения в формулу скорости

При чистом качении точка контакта в момент касания имеет нулевую скорость относительно земли. Значит, эта точка является мгновенным центром вращения твердого тела.

Если \(\theta(t)\) обозначает текущий угол ориентации эллипса, то центр движется вокруг точки контакта со скоростью

$$v(t)=\left|\frac{d\theta}{dt}\right|\lVert r(t)\rVert.$$

Тем самым геометрическая задача сводится к нахождению угловой скорости \(\theta'(t)\), возникающей из изменения направления касательной к эллипсу.

Шаг 3: Вычислить производную угла касательной

Для гладкой плоской параметризованной кривой производная угла касательной выражается формулой

$$\frac{d\theta}{dt}=\frac{x'(t)y''(t)-y'(t)x''(t)}{x'(t)^2+y'(t)^2}.$$

После подстановки производных эллипса получаем

$$x'(t)y''(t)-y'(t)x''(t)=ab,$$

$$x'(t)^2+y'(t)^2=a^2\sin^2 t+b^2\cos^2 t,$$

а значит

$$\frac{d\theta}{dt}=\frac{ab}{a^2\sin^2 t+b^2\cos^2 t}.$$

Для \(a,b>0\) знаменатель всегда положителен, поэтому угловая скорость также остается положительной, и в интеграле длины дуги не возникает неопределенности со знаком.

Шаг 4: Получить замкнутый интегранд скорости центра

Объединяя предыдущие шаги, получаем

$$v(t)=\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}.$$

Следовательно, полная длина пути центра за один оборот равна

$$C(a,b)=\int_0^{2\pi}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

Полезная проверка получается в круговом случае \(a=b=R\). Тогда интегранд вырождается в константу \(R\), и потому

$$C(R,R)=\int_0^{2\pi}R\,dt=2\pi R,$$

что полностью совпадает с ожидаемой длиной пути центра катящегося круга.

Шаг 5: Использовать симметрию и интегрировать только четверть оборота

Интегранд зависит только от \(\sin^2 t\) и \(\cos^2 t\). Поэтому

$$v(t+\pi)=v(t),\qquad v(\pi-t)=v(t).$$

Отсюда следует

$$C(a,b)=4\int_0^{\pi/2}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

Именно эту форму вычисляет реализация. С численной точки зрения она удобна, потому что интервал короткий, а сам интегранд гладок на \([0,\pi/2]\).

Шаг 6: Разобранный пример для \(C(2,4)\)

При \(a=2\) и \(b=4\) интегранд на четверти оборота принимает вид

$$v(t)=\frac{4\sqrt{\cos^2 t+4\sin^2 t}}{\sin^2 t+4\cos^2 t}.$$

В трех стандартных узлах формулы Симпсона получаем

$$v(0)=1,\qquad v\left(\frac{\pi}{4}\right)=\frac{4\sqrt{10}}{5},\qquad v\left(\frac{\pi}{2}\right)=8.$$

Одна грубая панель Симпсона на \([0,\pi/2]\) дает приближение

$$S=\frac{\pi/2}{6}\left(1+4\cdot\frac{4\sqrt{10}}{5}+8\right)=\frac{\pi}{12}\left(9+\frac{16\sqrt{10}}{5}\right)\approx 5.0056.$$

После умножения на \(4\) получается примерно \(20.0224\), чего недостаточно по точности. После адаптивного уточнения численный интеграл сходится к значению

$$C(2,4)\approx 21.38816906,$$

и это наглядно показывает, зачем нужен адаптивный метод вместо одной грубой панели.

Как работает код

Реализации на C++, Python и Java используют один и тот же алгоритм. Сначала вычисляется приведенная выше замкнутая формула скорости. Затем на четвертном интервале \([0,\pi/2]\) применяется правило Симпсона по значениям в концах и в середине.

Далее запускается адаптивная рекурсия. Для текущего интервала программа дополнительно вычисляет скорость в двух внутренних четвертных точках, сравнивает сумму левой и правой аппроксимаций Симпсона с аппроксимацией на целом интервале и трактует разность как индикатор ошибки.

Если локальная оценка уже удовлетворяет условию

$$|\Delta|\le 15\varepsilon,$$

или достигнут предел глубины рекурсии, интервал принимается и используется стандартная поправка Ричардсона

$$S_{\text{refined}}=S_{\text{left}}+S_{\text{right}}+\frac{\Delta}{15}.$$

Иначе интервал делится пополам, и тот же процесс повторяется на двух подинтервалах. Искомый интеграл получается умножением сошедшегося значения на четверти интервала на \(4\). После этого реализация вычисляет \(C(1,4)\) и \(C(3,4)\), складывает их и выводит результат с восемью знаками после запятой. Численные параметры таковы: \(\varepsilon=10^{-13}\) и максимальная глубина рекурсии \(30\).

Анализ сложности

У адаптивного метода Симпсона нет заранее фиксированного числа шагов, поэтому естественно выражать сложность через число принятых подинтервалов. Если рекурсия заканчивается на \(N\) принятых листьях, то суммарное время равно \(O(N)\) вычислений функции с точностью до постоянных множителей, потому что каждое разбиение добавляет лишь константное число новых точек. Память составляет \(O(D)\), где \(D\) — глубина рекурсии; в реализации выполнено \(D\le 30\).

В данной задаче интегранд гладкий на \([0,\pi/2]\), поэтому на практике \(N\) остается небольшим. Адаптивная стратегия эффективна именно потому, что она уточняет сетку только там, где этого требует форма интегранда, а не тратит одинаковое разрешение на весь интервал.

Примечания и ссылки

  1. Страница задачи: https://projecteuler.net/problem=525
  2. Ellipse: Wikipedia — Ellipse
  3. Длина дуги: Wikipedia — Arc length
  4. Адаптивный метод Симпсона: Wikipedia — Adaptive Simpson's method

ملخص المسألة

لدينا إهليلج أنصاف محاوره \(a\) و\(b\) يتدحرج على خط مستقيم دون انزلاق. نرمز بـ \(C(a,b)\) إلى طول المسار الذي يقطعه مركزه خلال دورة كاملة. المطلوب هو حساب

$$C(1,4)+C(3,4).$$

التنفيذ لا يبحث عن دالة أصلية ابتدائية مغلقة. بدلًا من ذلك يشتق أولًا صيغة دقيقة لسرعة المركز، ثم يحول المسألة الهندسية إلى تكامل أحادي البعد ناعم، وبعدها يحسب هذا التكامل عدديًا بدقة عالية.

المنهج الرياضي

نعمل في جهاز الإحداثيات الخاص بالإهليلج، ونُمثّل نقطة الحد التي تلامس الأرض في تلك اللحظة بالبارامتر

$$r(t)=(x(t),y(t))=(a\cos t,b\sin t),\qquad 0\le t\le 2\pi.$$

وعندما يمر \(t\) مرة واحدة على المجال الكامل، فإن نقطة التماس تمر على كامل محيط الإهليلج مرة واحدة أيضًا، ولذلك يمكن تمثيل دورة التدحرج الكاملة بالتكامل على دورة كاملة للمتغير \(t\).

الخطوة 1: بارمترة الإهليلج وكتابة هندسته المحلية

المشتقتان الأولى والثانية هما

$$r'(t)=(-a\sin t,b\cos t),\qquad r''(t)=(-a\cos t,-b\sin t).$$

أما المسافة من مركز الإهليلج إلى نقطة التماس فهي

$$\lVert r(t)\rVert=\sqrt{a^2\cos^2 t+b^2\sin^2 t}.$$

هذه المسافة ليست ثابتة إلا إذا كان \(a=b\). ولهذا السبب تحديدًا لا يتحرك مركز الإهليلج المتدحرج ببساطة مركز الدائرة المتدحرجة.

الخطوة 2: تحويل شرط التدحرج دون انزلاق إلى صيغة للسرعة

في التدحرج النقي تكون سرعة نقطة التماس بالنسبة إلى الأرض مساوية للصفر لحظة التماس، ولذلك فإن نقطة التماس تمثل مركز الدوران اللحظي للجسم الصلب.

إذا كان \(\theta(t)\) هو زاوية توجيه الإهليلج في تلك اللحظة، فإن سرعة المركز حول نقطة التماس تساوي

$$v(t)=\left|\frac{d\theta}{dt}\right|\lVert r(t)\rVert.$$

وهكذا تختزل المسألة الهندسية إلى إيجاد السرعة الزاوية \(\theta'(t)\) الناتجة عن تغير اتجاه المماس على الإهليلج.

الخطوة 3: اشتقاق زاوية المماس

لأي منحنى مستوٍ أملس معلم بارامتريًا، فإن مشتقة زاوية المماس تعطى بالعلاقة

$$\frac{d\theta}{dt}=\frac{x'(t)y''(t)-y'(t)x''(t)}{x'(t)^2+y'(t)^2}.$$

وبالتعويض من مشتقات الإهليلج نحصل على

$$x'(t)y''(t)-y'(t)x''(t)=ab,$$

$$x'(t)^2+y'(t)^2=a^2\sin^2 t+b^2\cos^2 t,$$

ومن ثم

$$\frac{d\theta}{dt}=\frac{ab}{a^2\sin^2 t+b^2\cos^2 t}.$$

وبما أن المقام موجب دائمًا عندما \(a,b>0\)، فإن السرعة الزاوية تبقى موجبة، ولا تبقى أي مشكلة إشارة داخل تكامل طول القوس.

الخطوة 4: الحصول على المدمج المغلق لسرعة المركز

بجمع النتيجتين السابقتين نحصل على

$$v(t)=\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}.$$

وعليه فإن طول مسار المركز خلال دورة كاملة يساوي

$$C(a,b)=\int_0^{2\pi}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

هناك فحص معقول مهم في حالة الدائرة \(a=b=R\). عندها ينهار المدمج إلى الثابت \(R\)، وبالتالي

$$C(R,R)=\int_0^{2\pi}R\,dt=2\pi R,$$

وهو بالضبط طول المسار المتوقع لمركز دائرة تتدحرج.

الخطوة 5: استخدام التناظر والاكتفاء بربع دورة

بما أن المدمج يعتمد فقط على \(\sin^2 t\) و\(\cos^2 t\)، فإن

$$v(t+\pi)=v(t),\qquad v(\pi-t)=v(t).$$

ومن هنا نستنتج

$$C(a,b)=4\int_0^{\pi/2}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

هذه هي الصيغة الدقيقة التي يحسبها التنفيذ فعليًا، وهي مناسبة عدديًا لأن المجال قصير ولأن المدمج ناعم على \([0,\pi/2]\).

الخطوة 6: مثال محلول لـ \(C(2,4)\)

عندما \(a=2\) و\(b=4\)، يصبح مدمج ربع الدورة

$$v(t)=\frac{4\sqrt{\cos^2 t+4\sin^2 t}}{\sin^2 t+4\cos^2 t}.$$

وعند عقد Simpson القياسية الثلاث نحصل على

$$v(0)=1,\qquad v\left(\frac{\pi}{4}\right)=\frac{4\sqrt{10}}{5},\qquad v\left(\frac{\pi}{2}\right)=8.$$

ولو استخدمنا لوحة Simpson واحدة فقط على \([0,\pi/2]\)، فسنحصل على التقدير الخشن

$$S=\frac{\pi/2}{6}\left(1+4\cdot\frac{4\sqrt{10}}{5}+8\right)=\frac{\pi}{12}\left(9+\frac{16\sqrt{10}}{5}\right)\approx 5.0056.$$

وبعد الضرب في \(4\) نحصل على نحو \(20.0224\)، وهي دقة غير كافية. أما بعد التنقيح التكيفي فإن التكامل العددي يتقارب إلى

$$C(2,4)\approx 21.38816906,$$

وهذا يوضح مباشرة لماذا استُخدمت الطريقة التكيفية بدلًا من لوحة Simpson خشنة واحدة.

كيف يعمل الكود

التنفيذات في C++ وPython وJava تتبع الخوارزمية نفسها. تبدأ أولًا بحساب صيغة السرعة المغلقة المذكورة أعلاه. بعد ذلك تطبق قاعدة Simpson على مجال ربع الدورة \([0,\pi/2]\) باستخدام القيمتين عند الطرفين والقيمة عند المنتصف.

ثم تبدأ العودية التكيفية. ففي أي مجال حالي، يحسب التنفيذ سرعة المركز عند نقطتي الربع الداخليتين، ثم يقارن مجموع تقريبي Simpson على النصفين مع تقريب Simpson على المجال الكامل، ويعتبر الفرق إشارة للخطأ المحلي.

إذا كان تقدير الخطأ يحقق بالفعل

$$|\Delta|\le 15\varepsilon,$$

أو إذا بلغنا حد العمق العودي، فإن الكود يقبل المجال ويطبق تصحيح Richardson القياسي

$$S_{\text{refined}}=S_{\text{left}}+S_{\text{right}}+\frac{\Delta}{15}.$$

وإلا فإنه يقسم المجال إلى نصفين ويكرر العملية نفسها على المجالين الفرعيين. بعد ذلك يضرب قيمة ربع الدورة المتقاربة في \(4\) للحصول على \(C(a,b)\). وفي النهاية يحسب التنفيذ \(C(1,4)\) و\(C(3,4)\)، ثم يجمعهما ويطبع الناتج بدقة ثمانية منازل عشرية. المعاملات العددية المستخدمة هي \(\varepsilon=10^{-13}\) وعمق عودي أقصى مقداره \(30\).

تحليل التعقيد

لا تمتلك طريقة Simpson التكيفية عددًا ثابتًا من الخطوات معروفًا مسبقًا، لذلك يكون التحليل الأنسب بدلالة عدد المجالات الفرعية المقبولة في النهاية. فإذا انتهت العودية عند \(N\) أوراق مقبولة، فإن الزمن الكلي يكون \(O(N)\) من تقييمات الدالة حتى العوامل الثابتة، لأن كل عملية تقسيم تضيف عددًا ثابتًا فقط من العينات الجديدة. أما الذاكرة فهي \(O(D)\)، حيث \(D\) هو عمق العودية، ومع التنفيذ هنا لدينا \(D\le 30\).

في هذه المسألة يكون المدمج ناعمًا على \([0,\pi/2]\)، ولذلك يبقى \(N\) صغيرًا عمليًا. والميزة الحقيقية للاستراتيجية التكيفية أنها تزيد الدقة فقط في المواضع التي يفرضها شكل المدمج، بدلًا من فرض نفس الكثافة الحسابية على كامل المجال.

حواشٍ ومراجع

  1. صفحة المسألة: https://projecteuler.net/problem=525
  2. Ellipse: Wikipedia — Ellipse
  3. طول القوس: Wikipedia — Arc length
  4. طريقة Simpson التكيفية: Wikipedia — Adaptive Simpson's method