I want to solve the following sets of $n$ coupled equations. Initial values of $x_{n}(t)$ and $p_{n}(t)$ are specified.

The problem is, I have an 1D lattice where every particle is bound with neighbouring particles using harmonic plus a slight nonlinear force depending on relative displacement of particles with the neighbouring ones.

$x_{n}(t)$, $p_{n}(t)$ are position and momentum of $n$-th particle.

The Hamilton's equations are as follows.

$$ \frac{d x_{n}(t)}{d t}=\frac{p_{n}(t)}{m} $$

$$ \frac{d p_{n}(t)}{d t}=-k ((x_{n}(t)-x_{n-1}(t))+(x_{n}(t)-x_{n+1}(t)))-a ((x_{n}(t)-x_{n-1}(t))^3+(x_{n}(t)-x_{n+1}(t))^3) $$

To solve this in the Runge-Kutta method, the main problem lies in the fact that all the $n$ equations are coupled and depend on not one, but three other variables, which depend on the other three and so on, where all of them depends on time. what would be the method for solving this using Runge-Kutta?

I tried with the steps mentioned in https://www.myphysicslab.com/explain/runge-kutta-en.html

but the result that came out was wrong. This is the code I ran.

!f=dp/dt=Force,force on nth particle=f(k,a,x[n-1],x[n],x[n+1]),
!1d lattice(2m+1) particles, all particles connected with each other in harmonic potential, with slight x^3 nonlinearity.

function f(k,a,q,y,z) result(s)
end function f

!g=dx/dt=p/mass(mass=1),p=momentum, particles at lattice ends cant move

function g(q) result(s)
end function g



write(*,*) "serial no of middlemost particle,time of observation"
read(*,*) m,tf


            !define initial momentum

                do i=1,2*m+1
                end do
            !define initial position(gaussian tilt in middle zone)

                do i1=1,m-10
                end do  

                do i2=m+10,2*m+1
                end do  

                do i3=m-9,m+9
                end do  


    do while(t.le.tf)  !stop at time of observation

        do i4=2,2*m

!(k1=first runge kutta coefficient for ith particle moving under force field of (i-1) and (i+1)th particle,k2=second and so on)
!started with i5=2, not 1, as, end particles are fixed, same for 2m+1th particle. k(1),k(2m+1) zero for same reason 

            do i5=2,2*m

            end do

            do i6=2,2*m

            end do

            do i7=2,2*m

            end do

            do i8=2,2*m

            end do


            x(i4)=x(i4)+h*(l1+2*l2+2*l3+l4)/(6*mass)        !final runge kutta result
        end do !for all particles(i4) at time t
    end do !for time=tf=time of observation

    !write to file


    do i9=1,2*m+1
    write(1,*) i9,x(i9)
    end do



The explicit code is not necessary, only the scheme for solving would help.

  • $\begingroup$ You tried to implement the RK by yourself right? Is it necessary to implement it by yourself or using a library like DifferentialEquations.jl is fine? $\endgroup$ Commented Jan 21, 2020 at 15:11
  • $\begingroup$ It would help to do myself as I have just started the numerical methods course. $\endgroup$
    – user157588
    Commented Jan 21, 2020 at 15:14
  • 1
    $\begingroup$ OK, so it's a homework right? so you need to put some other contexts here to say what do you mean exactly that your code doesn't work? Does it show any error? Results are wrong? If results are wrong, in what sense? Have you checked your implementation with a known analytical solution? Without these more information it's really difficult to help you. $\endgroup$ Commented Jan 21, 2020 at 15:18
  • $\begingroup$ Simply add an 'implicit none' (always do that!) to find errors in your code (maybe you are not declaring a program?). I think that it is not fair to give you a solved code, since this is a homework. I however modified yours to add some readability, hope this helps you have a clearer view. I suggest you to try it and to edit your question, but I will leave this to you: gist.github.com/GabrieleLabanca/…. $\endgroup$
    – Gabriele
    Commented Jan 23, 2020 at 11:50

3 Answers 3


To give a short summary of the answer of @Gabriele, the main point is that in a numerical solution of a second order system $x''=f(x)$ as you have with a first order solver as RK4, you need to transform the problem into a first order system $$ x'=p,\\ p'=f(x). $$ This is also the point of departure of the question, it is only in the implementation that the lines get crossed the wrong way.

So for instance in the simplest method, the Euler method, one would compute $$ x^{j+1}=x^j+hp^j,\\ p^{j+1}=p^j+hf(x^j). $$

For RK4, you would need to implement accordingly $\newcommand{\hide}[1]{\color{teal}{#1}}$ \begin{align} X^0&=x^j,&\hide{P^0}&\hide{=p^j,}&F^0=f(X^0),\\ X^1&=X^0+\tfrac12hP^0,&\hide{P^1}&\hide{=P^0+\tfrac12hF^0,}&F^1=f(X^1),\\ X^2&=X^0+\tfrac12hP^1,&\hide{P^2}&\hide{=P^0+\tfrac12hF^1,}&F^2=f(X^2),\\ X^3&=X^0+hP^2,&\hide{P^3}&\hide{=P^0+hF^2},&F^3=f(X^3),\\ \end{align} and then the step \begin{align} x^{j+1}&\hide{=x^j+\tfrac16h(P^0+2P^1+2P^2+P^3)}\\ &=x^j+hp^j+\tfrac16h^2(F^0+F^1+F^2),\\ p^{j+1}&=p^j+\tfrac16h(F^0+2F^1+2F^2+F^3).\\ \end{align} (the colored parts are not necessary for an implementation, they just highlight the steps from the standard RK4 method)


Runge-Kutta methods solve equations of the form $$\dot y = \frac{\mathrm{d} y}{\mathrm{d} t} = F(t,y(t)),$$ where $y$ can be multidimensional. The first step is reconducting your equation to this form. Starting from $$ \begin{cases} \dot x_n(t) = \frac{p_n(t)}{m} = f(p_{n}(t))& \\ \dot p_n(t) = -k \left[\left(x_n(t)-x_{n-1}(t)\right) + \left(x_n(t)-x_{n+1}(t)\right) \right] -a \left[\left(x_n(t)-x_{n-1}(t)\right)^3 + \left(x_n(t)-x_{n+1}(t)\right)^3 \right] = g(x_{\{n\}}(t))& \end{cases}$$ (where $x_{\{n\}}$ indicates that the function depends on many particles' positions) we can identify the multidimensional function $y = \{y_n\}_{n=1,\dots 2N}$ with $(\{x_n\},\{p_n\})_{n=1,\dots N}$ in the following way (note that the index refers to the particles, rather than to spacial dimension, as it would be with a vector): $$ \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \\ y_{N+1} \\ y_{N+2} \\ \vdots \\ y_{2N} \end{pmatrix} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_N \\ p_{1} \\ p_{2} \\ \vdots \\ p_{N} \end{pmatrix} $$ (note that $n$ runs to $2N$ for $y_n$, since it includes both the variables $x$ and $p$).

The time is discretized, so $t_i = t_0+i\cdot h$, and the corrisponding variables will be $y_n^i = y_n(t_i)$.

The equations finally become: $$ \begin{pmatrix} \dot y^i_n \\ \dot y^i_{2N+n} \\ \end{pmatrix} = \begin{pmatrix} \dot x^i_n \\ \dot p^i_{n} \\ \end{pmatrix} = \begin{pmatrix} \dfrac{p^i_n}{m} \\ -k \left[\left(x_n^i-x_{n-1}^i\right) + \left(x_n^i-x_{n+1}^i\right) \right] -a \left[\left(x^i_n-x^i_{n-1}\right)^3 + \left(x^i_n-x_{n+1}^i\right)^3 \right] \\ \end{pmatrix} = \begin{pmatrix} f_n\left(p_{\{n\}}^i\right) \\ g_n\left(x_{\{n\}}^i\right) \\ \end{pmatrix} = \begin{pmatrix} f_n^i \\ g_n^i \\ \end{pmatrix} $$

Now, as an example I will consider the RK2 (second order) method (see my thesis, Appendix A for the derivation); it is then easy to implement higher orders. The general scheme is $$ u_{i+1} = u_i + \frac{h}{2} F\left(t_i+\frac{h}{2},\, u_i+\frac{h}{2}\,F(t_i,\,u_i)\right) $$ and you must consider that: we are dealing with a multidimensional variable; the function has no explicit dependence on time; the form of the function, as well as the variables on which it depends, varies considering the positions or the momentums. Thus, $$ \begin{pmatrix} \dot x^{i+1}_{n} \\ \dot p^{i+1}_{n} \end{pmatrix} = \begin{pmatrix} \dot x^i_{n} + \frac{h}{2} f\left(p^i_n+\frac{h}{2}\cdot g_n^i\right)\\ \dot p^i_{n} + \frac{h}{2} g\left(x^i_n+\frac{h}{2}\cdot f_n^i\right) \end{pmatrix} $$ For example, for $x$, $$ x^{i+1}_n = x^i_n + \frac{h}{2} \frac{p^i_n+\frac{h}{2}\cdot g_n^i}{m} \\ = x^i_n + \frac{h}{2} \frac{p^i_n+\frac{h}{2}\cdot\left(-k \left[\left(x_n^i-x_{n-1}^i\right) + \left(x_n^i-x_{n+1}^i\right) \right] -a \left[\left(x^i_n-x^i_{n-1}\right)^3 + \left(x^i_n-x_{n+1}^i\right)^3 \right] \right) }{m} $$

To implement this, these steps are usual:

  • fill two arrays x0 and p0 with all the values at time $t_0$;
  • fill two arrays f0 and g0 evaluating the functions;
  • calculate f1 = f(p0+h/2.*g0) and g1 using all the arrays above;
  • calculate x1 = x0 + h/2.*f1 and p1.

You can then save, if necessary, the values and then assign the "1" values to the "0" arrays, in order to repeat the process.


I'm not sure why you think your system of ODEs cannot be solved by using RK4, but I think the easiest way to do this is using DifferentialEquations.jl. Basically you have $2n$ unknowns and $2n$ ODEs, which should be good to solve as long as you have their initial conditions. I think you could look at an example for system of ODEs here from their official website to find how you can solve it.

  • $\begingroup$ I have given the code by editing. $\endgroup$
    – user157588
    Commented Jan 21, 2020 at 15:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.