Computational Science Asked on August 22, 2021
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.
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 ways 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 begin{align} X^0&=x^j,&P^0&=p^j,&F^0=f(X^0), X^1&=X^0+tfrac12hP^0,&P^1&=P^0+tfrac12hF^0,&F^1=f(X^1), X^2&=X^0+tfrac12hP^1,&P^2&=P^0+tfrac12hF^1,&F^2=f(X^2), X^3&=X^0+hP^2,&P^3&=P^0+hF^2,&F^3=f(X^3), end{align} and then the step begin{align} x^{j+1}&=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}
Correct answer by Lutz Lehmann on August 22, 2021
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.
Answered by Alone Programmer on August 22, 2021
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+icdot 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}^iright) + left(x_n^i-x_{n+1}^iright) right] -a left[left(x^i_n-x^i_{n-1}right)^3 + left(x^i_n-x_{n+1}^iright)^3 right] end{pmatrix} = begin{pmatrix} f_nleft(p_{{n}}^iright) g_nleft(x_{{n}}^iright) 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} Fleft(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} fleft(p^i_n+frac{h}{2}cdot g_n^iright) dot p^i_{n} + frac{h}{2} gleft(x^i_n+frac{h}{2}cdot f_n^iright) 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}cdotleft(-k left[left(x_n^i-x_{n-1}^iright) + left(x_n^i-x_{n+1}^iright) right] -a left[left(x^i_n-x^i_{n-1}right)^3 + left(x^i_n-x_{n+1}^iright)^3 right] right) }{m} $$
To implement this, these steps are usual:
x0 and p0 with all the values at time $t_0$;f0 and g0 evaluating the functions;f1 = f(p0+h/2.*g0) and g1 using all the arrays above;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.
Answered by Gabriele on August 22, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP