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.