Aplicando o método Runge-Kutta aos ODEs de segunda ordem

11

Como posso substituir o método de Euler pela 4ª ordem de Runge-Kutta para determinar o movimento de queda livre em magnitude gravimétrica não constante (por exemplo, queda livre de 10.000 km acima do solo)?

Até agora escrevi uma integração simples pelo método Euler:

while()
{
    v += getMagnitude(x) * dt;
    x += v * dt;
    time += dt;
}

variável x significa posição atual, v significa velocidade, getMagnitude (x) retorna a aceleração na posição x.

Eu tentei implementar o RK4:

while()
{
    v += rk4(x, dt) * dt; // rk4() instead of getMagintude()
    x += v * dt;
    time += dt;
}

onde o corpo da função rk4 () é:

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;
}

Mas algo está errado, porque estou integrando apenas uma vez usando o RK4 (aceleração). A integração da velocidade usando o RK4 não faz sentido porque é o mesmo que v * dt.

Você poderia me dizer como resolver equações diferenciais de segunda ordem usando a integração Runge-Kutta? Devo implementar o RK4 calculando os coeficientes k1, l1, k2, l2 ... l4? Como eu posso fazer isso?

Marcin W.
fonte
Oi @ Marcin, editei o seu título para refletir melhor o que acho que realmente é o seu problema. Acho que podemos obter respostas mais úteis e será mais pesquisável para outras pessoas que virem essa pergunta no futuro com o novo título. Sinta-se à vontade para alterá-lo novamente se não concordar.
Doug Lipinski

Respostas:

17

Parece haver muita confusão sobre como aplicar métodos de várias etapas (por exemplo, Runge-Kutta) a ODEs de segunda ou mais alta ordem ou sistemas de ODEs. O processo é muito simples quando você o entende, mas talvez não seja óbvio sem uma boa explicação. O método a seguir é o que eu acho mais simples.

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))

Infelizmente, o C ++ não oferece suporte nativo a operações vetoriais como essa; portanto, você precisa usar uma biblioteca de vetores, usar loops ou gravar manualmente as partes separadas. Em C ++, você pode usar std::valarraypara obter o mesmo efeito. Aqui está um exemplo simples de trabalho com aceleração constante.

#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
fonte
6
" Infelizmente, o C ++ não suporta operações vetoriais como essa ", acho que sim, na biblioteca padrão, mas não necessariamente fácil de usar com outras bibliotecas de álgebra linear: en.cppreference.com/w/cpp/numeric/valarray , acho bibliotecas de álgebra linear comuns como Eigen, também devem contar como "suporte".
Kirill
1
@ Kirill, obrigado pela dica. Ainda sou relativamente novo em C ++ e nunca usei valarray antes, apenas aprendi algo útil também! Edição para adicionar.
Doug Lipinski
1
Talvez este conselho também seja útil: 1) Use o formato clang para formatar automaticamente seu código, é realmente padrão e uniforme. 2) Use typedef std::valarray<double> Vectorpara os tipos mais usados. 3) Use em const int NDIM = 2vez de #definepara o tipo de segurança e correção. 4) Desde o C ++ 11, você pode substituir o corpo do RHS simplesmente por return {X[1], 1}. 5) É realmente incomum em C ++ (ao contrário C) para variáveis primeiramente declarar, em seguida, mais tarde inicializar eles, preferem declarar variáveis no mesmo lugar onde você inicializa-la ( double t = 0., etc.)
Kirill
1
@MarcinW. RHS()calcula o lado direito da equação diferencial. O vetor de estado X é (x, v), portanto dX / dt = (dx / dt, dv / dt) = (v, a). Para o seu problema (se a = G * M / x ^ 2), o RHS deve retornar { X[1], G*M/(X[0]*X[0]) }.
Doug Lipinski
1
@ Kirill eu sei, mas isso só funciona desde o C ++ 11, o que significa que ele não funciona com as opções padrão do compilador nos compiladores mais populares. Eu escolhi deixar isso de fora em favor de algo que funcione com os padrões antigos também e espero reduzir a confusão causada pela incapacidade de compilar o código.
quer