1

Consider that I am solving a second order ODE using RK2/RK4. The ODE represents simple equations of motion: Equations of motion I am trying to solve:

\begin{align} \frac{dx}{dt} &= v \\[.3em] m·\frac{dv}{dt}&= f_{1}(x)+f_{2}(x,v) \end{align}

RK2 method:

\begin{align} s_{x1}&=h f_x (t_i, v_i) \\ s_{v1}&=h f_v (t_i, x_i, v_i) \\[.5em] s_{x2}&=h f_x(t_i+\tfrac{1}{2} h, v_i+\tfrac{1}{2}s_{v1}) \\ s_{v2}&=h f_v(t_i+\tfrac{1}{2}h, x_i+\tfrac{1}{2}s_{x1}, v_i +\tfrac{1}{2} s_{v1}) \\[1em] x_{i+1}&=x_i+s_{x2} \\ v_{i+1}&=v_{i}+s_{v2} \end{align} where $f_x(t,x,v)=v$ and $f_v(t,x,v)=\frac{1}{m}(f_1(x)+f_2(x,v))$

Now along with $x$ and $v$, I also require to compute the $f_{1}(x)+f_{2}(x,v)$ at each time-step $h$. In such case what should I take velocity and position pair at that particular time-step?

Method Runge Kutta $4^{th}$ order

Basic Formulae \begin{align} x^{'}&=x_{0}+ \frac{1}{6}(k_{0}+2k_{1}+2k_{2}+k_{3}) \\ v^{'}&=v_{0}+ \frac{1}{6}(l_{0}+2l_{1}+2l_{2}+l_{3}) \end{align}

Calculation of coefficients \begin{align} k_0 &= h v_0 \\ l_0 &= \frac {h (F_{p}(x_0) +F_g(x_0,v_0)) }{ m} \\[.5em] k_1 &= h (v_0+ \frac{l_0}{2}) \\ l_1 &= \frac {h (F_{p}(x_0 + \frac {k_0}{2}) +F_g(x_0 + \frac {k_0}{2}, v_0 + \frac {l_0}{2})) }{ m} \\[.5em] k_2 &= h (v_0+ \frac{ l_1}{2}) \\ l_2 &= \frac {h (F_{p}(x_0 + \frac {k_1}{2}) +F_g(x_0 + \frac {k_1}{2}, v_0 + \frac {l_1}{2})) }{ m} \\[.5em] k_3 &= h (v_0+ {l_2}) \\ l_3 &= \frac {h (F_{p}(x_0 + {k_2}) +F_g(x_0 + {k_2},v_0 + {l_2})) }{ m} \end{align}

In the above method, what will the value of $F_p$ and $F_g$. Should I take it the ones at the before applying runge-kutta simply at the $x_0$ and $v_0$ at that time-step. But this may seem incorrect as velocity and position are not computed using these force values.

Lutz Lehmann
  • 126,666
Abhishek Bhatia
  • 1,648

1 Answers1

1

You have to evaluate them exactly as the method prescribes. You already did this correctly for the RK2 method.

This might seem to be a lot more effort for RK4. But consider that RK4 is $O(h^4)$. Roughly, to get a accuracy of e.g. about $10^{-4}$ for $t=1$ you need $h=0.1$ and $10$ steps netting $40$ function evaluations. To get the same accuracy for the $O(h^2)$ RK2 method you need a step size $h=0.01$ and $100$ steps netting $200$ function evaluations.

A detailed example for using equal amounts of function evaluations in Euler, Heun, RK2, RK3 and RK4 (with one out of 3 in each method aiming for $10^{-4}$ accuracy) can be found in this answer: https://math.stackexchange.com/a/1239002/115115

Lutz Lehmann
  • 126,666
  • Thanks! But I am still a bit confused . Can you check the above edit please. – Abhishek Bhatia Apr 23 '15 at 11:42
  • Looks correct. I do not understand your final question. All the input arguments for the functions are perfectly defined, the functions can be evaluated without further preparation. – Lutz Lehmann Apr 23 '15 at 12:37
  • Using RK4 now at each time-step I can compute the velocity and position. But I tend to measure the forces too at each time-step. RK4 computes forces several times. Which value of force should I take? – Abhishek Bhatia Apr 23 '15 at 12:41
  • This is a different question. It depends. Is the function evaluation an experiment with freely selectable inputs or are you following a trajectory and trying to reconstruct it from its acceleration values measured in regular intervals? Runge-Kutta methods might be unusable in this situation, explore multi-step methods. – Lutz Lehmann Apr 23 '15 at 12:48
  • If I understand correctly the forces are computed directly and not using the acceleration values. So I wish to simplify measure the force at regular intervals of the timestep itself. – Abhishek Bhatia Apr 23 '15 at 13:05
  • But then you do not have a differential equation (only in the trivial sense). If you measure the acceleration $a(t)$ at regular intervals, $a_n=a(t_0+n·h)$, then the task is just to integrate twice using an appropriate integration scheme, Simpson or higher order. – Lutz Lehmann Apr 23 '15 at 13:15
  • Okay, I wish to measure individual forces $F_p$ and $F_g$. Acceleration will give me the total force. – Abhishek Bhatia Apr 23 '15 at 13:20
  • I think you should formulate exactly what data you have available, how the composition of that data set might be controlled (initial location, sampling frequency etc.) and how not, and what you wish to compute. As I understand, you wish to solve some kind of reverse problem, reconstruct or at least sample the functions $F_p$ and $F_g$ from some set of measured trajectories. – Lutz Lehmann Apr 23 '15 at 13:37
  • Okay, so I am writing simulation code using these equations. Consider sampling frequency to same as $h$. – Abhishek Bhatia Apr 23 '15 at 13:46