Solving coupled ODEs using Runge-Kutta method

Question

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)
real::k,a,q,y,z,s
s=-k*((y-q)+(y-z))-a*((y-q)**3+(y-z)**3)
end function f

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

function g(q) result(s)
real::s,q
s=q
end function g


real::x(1000),p(1000),k1(1000),k2(1000),k3(1000),k4(1000),h,mass,l1,l2,l3,l4,t,tf,a,k
integer::i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,m,i


mass=1
k=100
a=0
h=0.001




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

k=1000
a=10


            !define initial momentum

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

                do i1=1,m-10
                    x(i1)=0 
                end do  

                do i2=m+10,2*m+1
                    x(i2)=0 
                end do  

                do i3=m-9,m+9
                    x(i3)=exp(-(0.1*(i3-m))**2)
                end do  


    t=0

    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

                k1(1)=0
                k1(2*m+1)=0 
                k1(i5)=f(k,a,x(i5-1),x(i5),x(i5+1)) 
            end do


            do i6=2,2*m

                k2(1)=0
                k2(2*m+1)=0
                k2(i6)=f(k,a,x(i6-1)+1*h*k1(i6-1)/2,x(i6)+h*k1(i6)/2,x(i6+1)+1*h*k1(i6+1)/2)    
            end do


            do i7=2,2*m

                k3(1)=0
                k3(2*m+1)=0
                k3(i7)=f(k,a,x(i7-1)+1*h*k2(i7-1)/2,x(i7)+h*k2(i7)/2,x(i7+1)+1*h*k2(i7+1)/2)    
            end do


            do i8=2,2*m

                k4(1)=0
                k4(2*m+1)=0
                k4(i8)=f(k,a,x(i8-1)+1*h*k3(i8-1)/2,x(i8)+h*k3(i8)/2,x(i8+1)+1*h*k3(i8+1)/2)
            end do



            l1=g(p(i4))
            l2=g((p(i4)+h*l1/2))
            l3=g(p(i4)+h*l2/2)
            l4=g(p(i4)+h*l3)


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



    !write to file

        open(1,file="f.dat")

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

        close(1)

end

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

Answer 1

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)

Answer 2

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.

Answer 3

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.