Uji integrator symplectic orde-3 vs ke-4 dengan hasil yang aneh

Dalam jawaban saya untuk pertanyaan tentang MSE tentang simulasi fisika Hamiltonian 2D, saya telah menyarankan menggunakan integrator symplectic tingkat tinggi .

Kemudian saya pikir itu mungkin ide yang baik untuk menunjukkan efek dari langkah waktu yang berbeda pada keakuratan global metode dengan pesanan yang berbeda, dan saya menulis dan menjalankan skrip Python / Pylab untuk efek itu. Sebagai perbandingan, saya memilih:

• ( leap2 ) Contoh ke-2 Wikipedia yang saya kenal, meskipun saya tahu itu dengan nama leapfrog ,
• ( ruth3 ) integrator symplectic orde ketiga Ruth ,
• ( ruth4 ) Integrator symplectic orde 4 dari Ruth .

Yang aneh adalah, apa pun catatan waktu yang saya pilih, metode urutan ketiga Ruth tampaknya lebih akurat dalam pengujian saya daripada metode urutan keempat Ruth, bahkan dengan urutan besarnya.

Pertanyaan saya adalah: Apa yang saya lakukan salah di sini? Detail di bawah.

Metode membuka kekuatan mereka dalam sistem dengan Hamiltonians yang dapat dipisahkan , yaitu yang dapat ditulis sebagai mana terdiri dari semua koordinat posisi, terdiri dari momentum konjugat, mewakili kinetika konjugat, mewakili kinetik energi dan energi potensial

$$H(q,p) = T(p) + V(q)$$ $q$$p$$T$$V$

Dalam pengaturan kami, kami dapat menormalkan kekuatan dan momentum oleh massa yang digunakan. Dengan demikian kekuatan berubah menjadi akselerasi, dan momen berubah menjadi kecepatan.

Integrator symplectic datang dengan koefisien khusus (diberikan, konstan) yang akan saya label dan . Dengan koefisien tersebut, satu langkah untuk mengembangkan sistem dari waktu ke waktu berbentuk$a_1,\ldots,a_n$$b_1,\ldots,b_n$$t$$t+\delta t$

• Untuk :$i=1,\ldots,n$

  1. Hitung vektor dari semua akselerasi, diberi vektor dari semua posisi$g$$q$
  2. Ubah vektor dari semua kecepatan dengan$v$$b_i\,g\,\delta t$
  3. Ubah vektor dari semua posisi dengan$q$$a_i\,v\,\delta t$

Kebijaksanaan sekarang terletak pada koefisien. Ini adalah

$$\begin{align} \begin{bmatrix} a_1 & a_2 \\ b_1 & b_2 \end{bmatrix} &= \begin{bmatrix} \frac{1}{2} & \frac{1}{2} \\ 0 & 1 \end{bmatrix} &&\textsf{(leap2)} \\ \begin{bmatrix} a_1 & a_2 & a_3 \\ b_1 & b_2 & b_3 \end{bmatrix} &= \begin{bmatrix} \frac{2}{3} & -\frac{2}{3} & 1 \\ \frac{7}{24} & \frac{3}{4} & -\frac{1}{24} \end{bmatrix} &&\textsf{(ruth3)} \\ \begin{bmatrix} a_1 & a_2 & a_3 & a_4 \\ b_1 & b_2 & b_3 & b_4 \end{bmatrix} &= \frac{1}{2-\sqrt[3]{2}} \begin{bmatrix} \frac{1}{2} & \frac{1-\sqrt[3]{2}}{2} & \frac{1-\sqrt[3]{2}}{2} & \frac{1}{2} \\ 0 & 1 & -\sqrt[3]{2} & 1 \end{bmatrix} &&\textsf{(ruth4)} \end{align}$$

Untuk pengujian, saya telah memilih masalah nilai awal 1D yang memiliki Hamiltonian yang dapat dipisahkan. Di sini, diidentifikasi dengan .

$$\begin{align} y''+y &= 0 & y(0) &= 1 & y'(0) &= 0 \end{align}$$ $$\therefore\quad \left(y(t),y'(t)\right) = (\cos t, -\sin t)$$ $(q,v)$$(y,y')$

Saya telah mengintegrasikan IVP dengan metode di atas lebih dari dengan stepsize dari dengan integer dipilih di suatu tempat antara dan . Memperhatikan kecepatan leap2 , saya membuat tiga kali lipat untuk metode itu. Saya kemudian memplot kurva yang dihasilkan dalam ruang fase dan diperbesar di mana kurva idealnya harus kembali lagi pada .$t\in[0,2\pi]$$\delta t=\frac{2\pi}{N}$$N$$10$$100$$N$$(1,0)$$t=2\pi$

Berikut adalah plot dan zoom untuk dan :$N=12$$N=36$

Untuk , leap2 dengan ukuran langkah kebetulan lebih dekat ke rumah daripada ruth4 dengan ukuran langkah . Untuk , ruth4 menang atas leap2 . Namun, ruth3 , dengan ukuran langkah yang sama dengan ruth4 , tiba jauh lebih dekat ke rumah daripada yang lain, dalam semua pengaturan yang telah saya uji sejauh ini.$N=12$$\frac{2\pi}{3N}$$\frac{2\pi}{N}$$N=36$

Berikut ini skrip Python / Pylab:

import numpy as np
import matplotlib.pyplot as plt

def symplectic_integrate_step(qvt0, accel, dt, coeffs):
    q,v,t = qvt0
    for ai,bi in coeffs.T:
        v += bi * accel(q,v,t) * dt
        q += ai * v * dt
        t += ai * dt
    return q,v,t

def symplectic_integrate(qvt0, accel, t, coeffs):
    q = np.empty_like(t)
    v = np.empty_like(t)
    qvt = qvt0
    q[0] = qvt[0]
    v[0] = qvt[1]
    for i in xrange(1, len(t)):
        qvt = symplectic_integrate_step(qvt, accel, t[i]-t[i-1], coeffs)
        q[i] = qvt[0]
        v[i] = qvt[1]
    return q,v

c = np.math.pow(2.0, 1.0/3.0)
ruth4 = np.array([[0.5, 0.5*(1.0-c), 0.5*(1.0-c), 0.5],
                  [0.0,         1.0,          -c, 1.0]]) / (2.0 - c)
ruth3 = np.array([[2.0/3.0, -2.0/3.0, 1.0], [7.0/24.0, 0.75, -1.0/24.0]])
leap2 = np.array([[0.5, 0.5], [0.0, 1.0]])

accel = lambda q,v,t: -q
qvt0 = (1.0, 0.0, 0.0)
tmax = 2.0 * np.math.pi
N = 36

fig, ax = plt.subplots(1, figsize=(6, 6))
ax.axis([-1.3, 1.3, -1.3, 1.3])
ax.set_aspect('equal')
ax.set_title(r"Phase plot $(y(t),y'(t))$")
ax.grid(True)
t = np.linspace(0.0, tmax, 3*N+1)
q,v = symplectic_integrate(qvt0, accel, t, leap2)
ax.plot(q, v, label='leap2 (%d steps)' % (3*N), color='black')
t = np.linspace(0.0, tmax, N+1)
q,v = symplectic_integrate(qvt0, accel, t, ruth3)
ax.plot(q, v, label='ruth3 (%d steps)' % N, color='red')
q,v = symplectic_integrate(qvt0, accel, t, ruth4)
ax.plot(q, v, label='ruth4 (%d steps)' % N, color='blue')
ax.legend(loc='center')
fig.show()

Saya sudah memeriksa kesalahan sederhana:

• Tidak ada kesalahan ketik Wikipedia. Saya telah memeriksa referensi, khususnya ( 1 , 2 , 3 ).
• Saya sudah mendapatkan urutan koefisien yang benar. Jika Anda membandingkan dengan pemesanan Wikipedia, perhatikan bahwa urutan aplikasi operator berfungsi dari kanan ke kiri. Penomoran saya setuju dengan Candy / Rozmus . Dan jika saya mencoba pemesanan lain, hasilnya menjadi lebih buruk.

Kecurigaan saya:

• Urutan stepsize yang salah: Mungkin skema urutan ke-3 Ruth memiliki konstanta tersirat yang jauh lebih kecil, dan jika ukuran langkah dibuat sangat kecil, maka metode urutan ke-4 akan menang? Tetapi saya bahkan sudah mencoba , dan metode urutan ke-3 masih lebih unggul.$N=360$
• Tes salah: Sesuatu yang istimewa tentang pengujian saya membuat metode tingkat ketiga Ruth berperilaku seperti metode tingkat tinggi?

Mengikuti saran Kirill , saya menjalankan tes dengan dari daftar nilai yang secara geometris meningkat, dan untuk setiap menghitung kesalahan sebagai mana mewakili perkiraan diperoleh dengan integrasi numerik. Ini adalah hasil dari plot log-log:$N$$N$

$$\epsilon(N) = \|\tilde z(2\pi) - \tilde z(0)\|_2 \quad\text{where}\quad \tilde z(t) = \left(\tilde{y}(t),\tilde y'(t)\right)$$ $\tilde z$

Jadi ruth3 memang memiliki sama agar sebagai ruth4 untuk kasus uji dan konstanta tersirat hanya besarnya.$4$$\frac{1}{100}$

Menarik. Saya harus menyelidiki lebih lanjut, mungkin mencoba tes lain.

Ngomong-ngomong, berikut adalah tambahan skrip Python untuk plot kesalahan:

def int_error(qvt0, accel, qvt1, Ns, coeffs):
    e = np.empty((len(Ns),))
    for i,N in enumerate(Ns):
        t = np.linspace(qvt0[2], qvt1[2], N+1)
        q,v = symplectic_integrate(qvt0, accel, t, coeffs)
        e[i] = np.math.sqrt((q[-1]-qvt1[0])**2+(v[-1]-qvt1[1])**2)
    return e

qvt1 = (1.0, 0.0, tmax)
Ns = [12,16,20,24,32,40,48,64,80,96,128,160,192,
      256,320,384,512,640,768,1024,1280,1536,2048,2560,3072]

fig, ax = plt.subplots(1)
ax.set_xscale('log')
ax.set_xlabel(r"$N$")
ax.set_yscale('log')
ax.set_ylabel(r"$\|z(2\pi)-z(0)\|$")
ax.set_title(r"Error after 1 period vs #steps")
ax.grid(True)
e = int_error(qvt0, accel, qvt1, Ns, leap2)
ax.plot(Ns, e, label='leap2', color='black')
e = int_error(qvt0, accel, qvt1, Ns, ruth3)
ax.plot(Ns, e, label='ruth3', color='red')
e = int_error(qvt0, accel, qvt1, Ns, ruth4)
ax.plot(Ns, e, label='ruth4', color='blue')
ax.legend(loc='upper right')
fig.show()

Merencanakan kesalahan dari , selama interval penuh, diskalakan oleh kekuatan ukuran langkah yang diberikan oleh urutan yang diharapkan, memberikan plot$\ddot q=-q$$q(0)=1,\dot q(0)=0$

Seperti yang diharapkan, grafik untuk meningkatkan jumlah sub-interval semakin mendekati kurva batas yang merupakan koefisien kesalahan utama. Dalam semua kecuali satu plot konvergensi ini terlihat sangat cepat, hampir tidak ada perbedaan. Ini berarti bahwa bahkan untuk ukuran langkah yang relatif besar istilah kesalahan utama mendominasi semua istilah lainnya.

Dalam metode Rut urutan ke-3, koefisien utama komponen tampaknya nol, kurva batas terlihat dekat atau sama dengan sumbu horizontal. Grafik yang terlihat jelas menunjukkan dominasi istilah kesalahan urutan keempat. Penskalaan untuk kesalahan urutan ke-4 memberikan plot yang mirip dengan yang lain.$p$

Seperti yang dapat kita lihat, dalam semua 3 kasus koefisien kesalahan urutan utama dalam komponen adalah nol pada setelah periode penuh. Dalam menggabungkan kesalahan dari kedua komponen, perilaku komponen mendominasi, memberikan kesan metode urutan ke-4 dalam plot loglog.$q$$t=2\pi$$p$

Koefisien maksimal dalam komponen dapat ditemukan sekitar . Plot loglog di sana harus mencerminkan urutan kesalahan global yang benar.$q$$3\pi/2$

---

Bahwa lenyapnya istilah kesalahan tingkat 3 dalam Ruth3p adalah artefak dari kesederhanaan persamaan linear menunjukkan contoh non-linear , dengan plot yang sesuai$\ddot q=-\sin(q)$$q(0)=1.3,~\dot q(0)=0$