Menerapkan metode Runge-Kutta ke ODE urutan kedua

11

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?

Marcin W.
sumber
Hai @ Marsin, saya mengedit judul Anda untuk lebih mencerminkan apa yang menurut saya masalah Anda sebenarnya. Saya pikir kami mungkin mendapatkan jawaban yang lebih bermanfaat dan akan lebih mudah dicari orang lain yang melihat pertanyaan ini di masa depan dengan judul yang baru. Jangan ragu untuk mengubahnya kembali jika Anda tidak setuju.
Doug Lipinski

Jawaban:

17

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.

F=mx¨

[x˙v˙]=[vF/m]

vxk1k4(x,v)

while (t<TMAX)
    k1 = RHS( t, X );
    k2 = RHS( t + dt / 2, X + dt / 2 * k1 );
    k3 = RHS( t + dt / 2, X + dt / 2 * k2 );
    k4 = RHS( t + dt, X + dt * k3 );
    X = X + dt / 6 * ( k1 + 2 * k2 + 2 * k3 + k4 );
    t = t + dt;
end

X=(x,v)RHS( t, X )(x˙(t),v˙(t))

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 menggunakan std::valarrayuntuk mencapai efek yang sama. Inilah contoh kerja sederhana dengan akselerasi konstan.

#include <valarray>
#include <iostream>

const size_t NDIM = 2;

typedef std::valarray<double> Vector;

Vector RHS( const double t, const Vector X )
{
  // Right hand side of the ODE to solve, in this case:
  // d/dt(x) = v;
  // d/dt(v) = 1;
  Vector output(NDIM);
  output[0] = X[1];
  output[1] = 1;
  return output;
}

int main()
{

  //initialize values

  // State variable X is [position, velocity]
  double init[] = { 0., 0. };
  Vector X( init, NDIM );

  double t = 0.;
  double tMax=5.;
  double dt = 0.1;

  //time loop
  int nSteps = round( ( tMax - t ) / dt );
  for (int stepNumber = 1; stepNumber<=nSteps; ++stepNumber)
  {

    Vector k1 = RHS( t, X );
    Vector k2 = RHS( t + dt / 2.0,  X + dt / 2.0 * k1 );
    Vector k3 = RHS( t + dt / 2.0, X + dt / 2.0 * k2 );
    Vector k4 = RHS( t + dt, X + dt * k3 );

    X += dt / 6.0 * ( k1 + 2.0 * k2 + 2.0 * k3 + k4 );
    t += dt;
  }
  std::cout<<"Final time: "<<t<<std::endl;
  std::cout<<"Final position: "<<X[0]<<std::endl;
  std::cout<<"Final velocity: "<<X[1]<<std::endl;

}
Doug Lipinski
sumber
6
" Sayangnya C ++ tidak secara native mendukung operasi vektor seperti ini " Saya pikir memang demikian, bahkan di perpustakaan standar, tetapi tidak selalu mudah digunakan dengan perpustakaan aljabar linier lainnya: en.cppreference.com/w/cpp/numeric/valarray Saya pikir perpustakaan aljabar linear umum seperti Eigen, juga harus dihitung sebagai "dukungan".
Kirill
1
@ Kirill, Terima kasih atas tipnya. Saya masih relatif baru untuk C ++ dan saya belum pernah menggunakan valarray sebelumnya, saya baru belajar sesuatu yang berguna juga! Editing untuk ditambahkan.
Doug Lipinski
1
Mungkin saran ini juga akan sangat membantu kemudian: 1) Gunakan format dentang untuk memformat kode Anda secara otomatis, itu benar-benar standar dan seragam. 2) Gunakan typedef std::valarray<double> Vectoruntuk jenis yang umum digunakan. 3) Gunakan const int NDIM = 2sebagai ganti #defineuntuk keamanan dan kebenaran tipe. 4) Karena C ++ 11 Anda dapat mengganti tubuh RHS hanya dengan return {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.)
Kirill
1
@MarcinW. 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]) }.
Doug Lipinski
1
@ Kirill saya tahu, tapi itu hanya bekerja sejak C ++ 11 yang berarti tidak bekerja dengan opsi kompiler default pada kompiler paling populer. Saya memilih untuk mengabaikannya demi sesuatu yang bekerja dengan standar lama juga dan mudah-mudahan mengurangi kebingungan yang disebabkan oleh ketidakmampuan untuk menyusun kode.
Doug Lipinski