Bagaimana saya bisa mengganti metode Euler dengan urutan Runge-Kutta ke-4 untuk menentukan gerakan jatuh bebas dalam besaran gravitasi tidak konstan (mis. Jatuh bebas dari 10.000 km di atas tanah)?
Sejauh ini saya menulis integrasi sederhana dengan metode Euler:
while()
{
v += getMagnitude(x) * dt;
x += v * dt;
time += dt;
}
x variabel berarti posisi saat ini, v berarti kecepatan, getMagnitude (x) mengembalikan akselerasi pada posisi x.
Saya mencoba menerapkan RK4:
while()
{
v += rk4(x, dt) * dt; // rk4() instead of getMagintude()
x += v * dt;
time += dt;
}
di mana fungsi rk4 () adalah:
inline double rk4(double tx, double tdt)
{
double k1 = getMagnitude(tx);
double k2 = getMagnitude(tx + 0.5 * tdt * k1);
double k3 = getMagnitude(tx + 0.5 * tdt * k2);
double k4 = getMagnitude(tx + tdt * k3);
return (k1 + 2*k2 + 2*k3 + k4)/6.0;
}
Tapi ada yang salah, karena saya hanya mengintegrasikan sekali menggunakan RK4 (akselerasi). Mengintegrasikan kecepatan menggunakan RK4 tidak masuk akal karena sama dengan v * dt.
Bisakah Anda memberi tahu saya cara menyelesaikan persamaan diferensial orde kedua menggunakan integrasi Runge-Kutta? Haruskah saya menerapkan RK4 dengan menghitung k1, l1, k2, l2 ... l4 koefisien? Bagaimana saya bisa melakukan itu?
sumber
Jawaban:
Tampaknya ada sedikit kebingungan tentang bagaimana menerapkan metode multi-langkah (misalnya Runge-Kutta) ke ODE ke-2 atau lebih tinggi atau sistem ODE. Prosesnya sangat sederhana begitu Anda memahaminya, tetapi mungkin tidak jelas tanpa penjelasan yang baik. Metode berikut adalah yang saya temukan paling sederhana.
k1
k4
X
RHS( t, X )
Sayangnya C ++ tidak mendukung operasi vektor seperti ini sehingga Anda perlu menggunakan perpustakaan vektor, menggunakan loop, atau secara manual menuliskan bagian-bagian yang terpisah.Di C ++ Anda dapat menggunakanstd::valarray
untuk mencapai efek yang sama. Inilah contoh kerja sederhana dengan akselerasi konstan.sumber
typedef std::valarray<double> Vector
untuk jenis yang umum digunakan. 3) Gunakanconst int NDIM = 2
sebagai ganti#define
untuk keamanan dan kebenaran tipe. 4) Karena C ++ 11 Anda dapat mengganti tubuh RHS hanya denganreturn {X[1], 1}
. 5) Benar-benar tidak biasa dalam C ++ (tidak seperti C) untuk pertama-tama mendeklarasikan variabel, kemudian menginisialisasi mereka, lebih suka mendeklarasikan variabel di tempat yang sama di mana Anda menginisialisasi mereka (double t = 0.
, dll.)RHS()
menghitung Sisi Kanan dari persamaan diferensial. Vektor state X adalah (x, v) jadi dX / dt = (dx / dt, dv / dt) = (v, a). Untuk masalah Anda (jika a = G * M / x ^ 2) RHS harus kembali{ X[1], G*M/(X[0]*X[0]) }
.