Modelling an ODE

New to Julia
,
1

Dear Julia users,
I am trying to simulate the behaivour of the particles interacting with Lennard - Jones potential, but, at last, it seems to be not correct.

using DifferentialEquations, Plots, Interact
t=t+1
N=20 #Number of Particles
indexN = 4*N
end_time = 1 # Final time for the simulation

function normv(u1,u2) 
   return (u1^2+u2^2) 
end
function PotGradient(r) # The gradient of the interaction potential U(r) w.r.t r
    return 4*(12*(1/r)^13-6*(1/r)^7)
end
function Potential(u,index) # Total interaction potential
    # compute 
    pot =[0,0]                                                                
    coord =[u[index],u[index+1]] # =Xᵢ (position of particle 'i')         
    for j=1:4:indexN # run over all the particles
        if j!=index
            coordj = [u[j],u[j+1]] # =Xⱼ (position of particle 'j')
            pot+=[PotGradient(u[index]-u[j]), PotGradient(u[index+1]-u[j+1])]
        end
    end
    return pot  
end

function LJ(du,u,p,t)
    pot_aux =[0.0, 0.0]
    for i=1:4:indexN
        du[i]=u[i+2]  # dX_i/dt =V_i (first component)
        du[i+1] =u[i+3]  # dX_i/dt = V_i (second component)  # κ *∇ᵢU/N (Xᵢ)   
        pot_aux = Potential(u,i,p)  # κ *∇ᵢU/N (Xᵢ)  
        du[i+2]= pot_aux[1]
        du[i+3]= pot_aux[2]
    end
end

u_init=rand(4*N)*10
probCS = ODEProblem(LJ,u_init,(0.0,end_time))
solCS =solve(probCS)

retcode: Success
Interpolation: automatic order switching interpolation
t: 743-element Vector{Float64}:
 0.0
 2.1812646110003302e-12
 3.215320872725304e-12
 5.71276172412876e-12
 8.723722211746462e-12
 1.281302573566052e-11
 1.844953612690731e-11
 2.78878489735781e-11
 3.873403731143217e-11
 5.453440457066641e-11
 7.633119324798309e-11
 1.060160285219362e-10
 1.4469125904201344e-10
 ⋮
 1.803717165348969e-5
 1.8316562858629184e-5
 1.8881220311340966e-5
 1.9658047154716174e-5
 2.1346997254121142e-5
 2.5607873176011568e-5
 4.653112595381582e-5
 0.0002557636537318583
 0.0023480889315122826
 0.023271341709316527
 0.23250386948735896
 1.0
u: 743-element Vector{Vector{Float64}}:
 [5.344975321034342, 2.178050884532454, 4.259753183115283, 1.4431736260409322, 0.015634617932931416, 2.6694030821274906, 9.768506957862531, 3.653635041775982, 0.6447160727350965, 5.035146303684259  …  1.1806964139862797, 1.1694289721974749, 1.0233383910068894, 4.004042094274527, 9.537981839246727, 9.362649805515623, 0.4891350080035828, 4.3343658390368205, 3.461108251111551, 9.36348814355907]
 [5.344989187507778, 2.178050884535602, 1.2698431589318363e7, 1.4428830811383466, 0.01563461795423911, 2.6694030821364616, 9.768505183280789, 4.571753265805987, 0.6447160727465618, 5.035146303696906  …  1.1808504402829991, -8.49203469016822e6, 1.0233383910278453, 4.004042094295535, 9.676474857019098, 9.8996517234608, 0.4891350080074831, 4.33423158023816, 0.11508713674839921, -1.2137599229940319e8]
 [5.3450054073500795, 2.178050884537094, 1.866422728954793e7, 1.4427453446279068, 0.015634617964340295, 2.669403082141414, 9.768504342017719, 5.006998899567643, 0.6447160727541598, 5.035146303702901  …  1.1809234922891048, -1.2494313004832864e7, 1.0233383910378853, 4.004042094305904, 9.742129235420569, 10.15422391082523, 0.489135008006782, 4.334078781076774, -1.4711367744530002, -1.733214272706606e8]
 [5.345069753485448, 2.178050884540697, 3.2784485435694262e7, 1.4424126849838455, 0.015634617988736554, 2.6694030821552315, 9.768502310208778, 6.058199173158146, 0.6447160727782468, 5.035146303717381  …  1.1811001109493395, -2.20343167899132e7, 1.0233383910624136, 4.004042094332031, 9.900696938823554, 10.769061612996714, 0.48913500799832405, 4.333511526879466, -5.302166608693582, -2.749950892428887e8]
 [5.3451931279768585, 2.1780508845450393, 4.898837835254502e7, 1.4420116244166181, 0.015634618018149124, 2.6694030821753802, 9.768499860622656, 7.325545497373942, 0.6447160728180746, 5.035146303734838  …  1.181313591878053, -3.3170687957181834e7, 1.023338391092512, 4.004042094365572, 10.091869069661913, 11.510315369388993, 0.489135007975406, 4.332555627665377, -9.920926410761874, -3.525861943252876e8]
 [5.3454352481276075, 2.178050884550935, 6.900134439517802e7, 1.441466928330853, 0.01563461805809548, 2.669403082208856, 9.768496533743669, 9.046778238833802, 0.644716072891052, 5.035146303758547  …  1.1816049799081219, -4.737463306486175e7, 1.0233383911343117, 4.0040420944147, 10.351507437831424, 12.51702744924266, 0.48913500792201037, 4.330993089210022, -16.193845190669048, -4.041110969425047e8]
 [5.345891713996445, 2.178050884559058, 9.20535672147273e7, 1.4407161439653695, 0.0156346181131557, 2.6694030822665344, 9.768491948124641, 11.419247321837995, 0.6447160730272977, 5.035146303791227  …  1.1820105106701655, -6.471988218243942e7, 1.0233383911936667, 4.004042094489163, 10.709381180891256, 13.904602117327974, 0.4891350078063661, 4.328630199956869, -24.840151617026695, -4.2925219871603274e8]
 [5.346898149266555, 2.1780508845726496, 1.1894070693388698e8, 1.4394589587124365, 0.015634618205353747, 2.669403082393061, 9.76848426952489, 15.391936737329278, 0.6447160733479862, 5.0351463038459485  …  1.1827036157490178, -8.726294337311208e7, 1.023338391297573, 4.004042094631364, 11.308639251534382, 16.22799271200854, 0.48913500750359207, 4.324521038063816, -39.318354700769504, -4.386482738283898e8]
 [5.348291608692611, 2.1780508845882545, 1.3609550447920814e8, 1.438014244140811, 0.01563461831130452, 2.6694030825847626, 9.768475445538915, 19.957216296623876, 0.644716073859614, 5.035146303908832  …  1.1835285056901834, -1.0410123774794121e8, 1.0233383914239633, 4.0040420948218545, 11.997286258192918, 18.897773796457034, 0.4891350069869091, 4.319752033720185, -55.956212466572836, -4.402531226955847e8]
 [5.3505426762334265, 2.178050884610959, 1.4701516666502044e8, 1.4359096318727458, 0.015634618465649917, 2.6694030829526345, 9.768462591048644, 26.60776381862871, 0.6447160748787893, 5.035146304000439  …  1.18479809232532, -1.1699587391082464e8, 1.0233383916214502, 4.00404209515117, 13.000484365589493, 22.786579223843162, 0.4891350059112999, 4.312792683074943, -80.19368846301654, -4.4054837176853544e8]
 [5.353809788751962, 2.1780508846422255, 1.5172900235101384e8, 1.4330062949460778, 0.015634618678570838, 2.6694030836325857, 9.768444858130392, 35.78227008656195, 0.6447160768177893, 5.035146304126812  …  1.1867094452415567, -1.2388017857231775e8, 1.0233383919199015, 4.0040420957063, 14.384407715213479, 28.149996734733577, 0.4891350037989377, 4.303189612462294, -113.6295574144877, -4.405829017897967e8]
 [5.358341357224842, 2.1780508846847053, 1.531911120941129e8, 1.4290522683801548, 0.015634618968545157, 2.6694030848802277, 9.76842070783877, 48.276939570670876, 0.6447160804525155, 5.035146304298918  …  1.1896760852156094, -1.2653464456025046e8, 1.0233383923748745, 4.004042096650299, 16.269159525515214, 35.450471735625065, 0.48913499974999913, 4.290110918057826, -159.1655213636824, -4.40585866358387e8]
 [5.3642743263728745, 2.1780508847398745, 1.5351926201214704e8, 1.4239007188629906, 0.015634619346340473, 2.6694030870621432, 9.768389243347437, 64.55576003437423, 0.6447160869077249, 5.035146304523148  …  1.1943340222444319, -1.27251112386171e8, 1.023338393051573, 4.004042098205084, 18.724730229550328, 44.947792625660306, 0.4891349924469929, 4.273071152543529, -218.49252554650232, -4.405860640136315e8]
 ⋮
 [9.85423244923324, -3841.4790909577673, 428927.13127142156, -2.1321216918011922e8, -2138.3596839015822, -2769.426419615532, -1.1876827200209141e8, -1.5374626237385026e8, -139.92737684848422, -632.6915171170268  …  1.4072840008144083e6, -3.268792193362076e6, -133.27190231930666, -2107.9912862937213, -7.429739418945514e6, -1.1699171318745528e8, -145.9112272065992, -2039.5427675921949, -8.196666134229385e6, -1.1339913883508603e8]
 [9.974070917356187, -3901.048695855406, 428927.1312714164, -2.1321216918011922e8, -2171.5424945485815, -2812.3817731458544, -1.1876827200209141e8, -1.5374626237385026e8, -142.10321788838434, -642.5874940899602  …  1.4072840008144083e6, -3.268792193362076e6, -135.34770616943834, -2140.6777420324984, -7.429739418946156e6, -1.1699171318745528e8, -148.2013036359676, -2071.22548965311, -8.196666134229385e6, -1.1339913883508603e8]
 [10.216267818698874, -4021.4405361918057, 428927.1312714091, -2.1321216918011922e8, -2238.6058844802624, -2899.19574602183, -1.1876827200209141e8, -1.5374626237385026e8, -146.50065358605545, -662.587540787354  …  1.4072840008144083e6, -3.268792193362076e6, -139.54296390405315, -2206.7379847893144, -7.429739418947367e6, -1.1699171318745528e8, -152.8296122560505, -2135.2571585274395, -8.196666134229385e6, -1.1339913883508603e8]
 [10.549469928122424, -4187.0694725451785, 428927.1312714033, -2.1321216918011922e8, -2330.868266312775, -3018.6299697024447, -1.1876827200209141e8, -1.5374626237385026e8, -152.5504207069208, -690.1025821078543  …  1.4072840008144083e6, -3.268792193362076e6, -145.31458492397496, -2297.620288045783, -7.429739418948804e6, -1.1699171318745528e8, -159.19700253530436, -2223.3486535901666, -8.196666134229385e6, -1.1339913883508603e8]
 [11.273906449120743, -4547.17418687629, 428927.1312713979, -2.1321216918011922e8, -2531.461951116864, -3278.2997348219014, -1.1876827200209141e8, -1.5374626237385026e8, -165.70361548552506, -749.9248348655968  …  1.4072840008144083e6, -3.268792193362076e6, -157.86304405416362, -2495.2134536632934, -7.4297394189510625e6, -1.1699171318745528e8, -173.04076261750043, -2414.874140398123, -8.196666134229385e6, -1.1339913883508603e8]
 [13.101511735000567, -5455.644784789888, 428927.13127139525, -2.1321216918011922e8, -3037.5188215751077, -3933.3934822512874, -1.1876827200209141e8, -1.5374626237385026e8, -198.88643772984986, -900.8441524099754  …  1.4072840008144083e6, -3.268792193362076e6, -189.52024185029882, -2993.700627444432, -7.429739418953588e6, -1.1699171318745528e8, -207.96573998861308, -2898.05380062365, -8.196666134229385e6, -1.1339913883508603e8]
 [22.076062525850382, -9916.736895849486, 428927.1312713948, -2.1321216918011922e8, -5522.537398657877, -7150.265393541969, -1.1876827200209141e8, -1.5374626237385026e8, -361.8324118503346, -1641.941308938709  …  1.4072840008144083e6, -3.268792193362076e6, -344.97455778630246, -5441.547815373933, -7.429739418954845e6, -1.1699171318745528e8, -379.466657450362, -5270.732647255472, -8.196666134229385e6, -1.1339913883508603e8]
 [111.82157043434852, -54527.65800644546, 428927.1312713948, -2.1321216918011922e8, -30372.723169485565, -39318.98450644879, -1.1876827200209141e8, -1.5374626237385026e8, -1991.2921530551816, -9052.912874226044  …  1.4072840008144083e6, -3.268792193362076e6, -1899.517717146406, -29920.019694668943, -7.429739418954928e6, -1.1699171318745528e8, -2094.475832067851, -28997.521113573697, -8.196666134229385e6, -1.1339913883508603e8]
 [1009.2766495193299, -500636.8691124051, 428927.1312713948, -2.1321216918011922e8, -278874.5808777624, -361006.17563551693, -1.1876827200209141e8, -1.5374626237385026e8, -18285.88956510365, -83162.62852709938  …  1.4072840008144083e6, -3.268792193362076e6, -17444.94931074744, -274704.738487619, -7.429739418954928e6, -1.1699171318745528e8, -19244.56757824274, -266265.4057767559, -8.196666134229385e6, -1.1339913883508603e8]
 [9983.827440369145, -4.961728980172002e6, 428927.1312713948, -2.1321216918011922e8, -2.763893157960531e6, -3.577878086926198e6, -1.1876827200209141e8, -1.5374626237385026e8, -181231.86368558835, -824259.7850558327  …  1.4072840008144083e6, -3.268792193362076e6, -172899.26524675777, -2.7225519264171193e6, -7.429739418954928e6, -1.1699171318745528e8, -190745.48503999162, -2.638944252408578e6, -8.196666134229385e6, -1.1339913883508603e8]
 [99729.33534886729, -4.9572650090767965e7, 428927.1312713948, -2.1321216918011922e8, -2.7614078928788215e7, -3.574659719983301e7, -1.1876827200209141e8, -1.5374626237385026e8, -1.8106916048904352e6, -8.235231350343166e6  …  1.4072840008144083e6, -3.268792193362076e6, -1.727442424606861e6, -2.7201023805712122e7, -7.429739418954928e6, -1.1699171318745528e8, -1.9057546596574804e6, -2.63657327187268e7, -8.196666134229385e6, -1.1339913883508603e8]
 [428929.2488715504, -2.13212164914716e8, 428927.1312713948, -2.1321216918011922e8, -1.1876826811806619e8, -1.537462586525443e8, -1.1876827200209141e8, -1.5374626237385026e8, -7.78779305982058e6, -3.541978080878615e7  …  1.4072840008144083e6, -3.268792193362076e6, -7.429738679372005e6, -1.1699171097912875e8, -7.429739418954928e6, -1.1699171318745528e8, -8.196664200782539e6, -1.1339913297812106e8, -8.196666134229385e6, -1.1339913883508603e8]

anim = @animate for time ∈ 0.0:0.1:end_time
      s = solCS(time)
    x=[]
    y=[]
    u=[]
    v=[]
    for i=1:4:indexN
        append!(x,s[i])
        append!(y,s[i+1])
        append!(u,(s[i+2]))
        append!(v,(s[i+3]))
    end
    fig=scatter(x,y,xlims=[0,2],ylims=[0,2]) 
    quiver!(x,y,quiver=(0.1*u,0.1*v)) 
     
end
    gif(anim)

And then there is no gif.
How could I improve the code?
Thanks.

2

I think you want push! rather than append!

3

I had to make a few changes to get this to run, I’m not sure what this line does, so I deleted it:

t=t+1

When you call Potential, you call a 3-argument method that does not exist, so I changed

        pot_aux = Potential(u,i,p)  # κ *∇ᵢU/N (Xᵢ)

to

        pot_aux = Potential(u,i)  # κ *∇ᵢU/N (Xᵢ)

With the parameters you selected, the solution blows up very quickly

julia> maximum(abs.(solCS(0.01)[1:4:end]))
9.078799476865399e7

julia> maximum(abs.(solCS(0.01)[2:4:end]))
191577.2029134577

So I plotted only very early times and set some extreme plot limits

anim = @animate for time ∈ 0.0:0.01:0.1
    s = solCS(time)
    x=[]
    y=[]
    u=[]
    v=[]
    for i=1:4:indexN
        append!(x,s[i])
        append!(y,s[i+1])
        append!(u,s[i+2])
        append!(v,s[i+3])
    end
    fig=scatter(x,y,xlims=[-1e8,1e8],ylims=[-1e5,1e5])
    quiver!(x,y,quiver=(0.1*u,0.1*v))
end

It is also helpful to set a name for the file to save.

gif(anim, "LJ.gif")

LJ

4

Thank you for answer, but using your advices (with other parameters) I got the following error while computing the solution:

┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.

└ @ SciMLBase C:\Users\moses\.julia\packages\SciMLBase\nbKmA\src\integrator_interface.jl:345

Out[189]:

retcode: DtLessThanMin

Interpolation: automatic order switching interpolation

t: 1-element Vector{Float64}:

0.0

u: 1-element Vector{Vector{Float64}}:

[6.718121762315706, 6.65602508203267, 2.5920043182956176, 0.05046174816268456, 4.114756258255367, 5.788530455588443, 3.63170289121328, 9.592034575664304, 7.025454672154869, 1.40231594907122 … 0.4255816536847723, 9.178614618856173, 2.7388950341671348, 7.960459970988758, 7.519797966139256, 2.5887995078875514, 2.596483542155319, 7.263343824685418, 2.4114191829104947, 0.3697416414917498]

Is the reason the initial condition?
How could I fix it?
Update:
Sometimes I can get solution, but animation still doesn’t work

5

This is a common error with many possible causes, first give the FAQ a look. If that doesn’t help, post exactly the code that gives the error. Try setting a random seed to ensure the example can be reproduced.