by

Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Activity
  • Sep 18 09:08

    github-actions[bot] on v1.15.1

    (compare)

  • Sep 18 08:53

    apbraun on exp-fix

    (compare)

  • Sep 18 08:41
    JuliaRegistrator commented on 440a1f2
  • Sep 18 08:40
    Datseris commented on 440a1f2
  • Sep 18 08:40

    Datseris on master

    Changed `e` to `Base.MathConsta… (compare)

  • Sep 18 08:40
    Datseris closed #125
  • Sep 18 06:59
    apbraun synchronize #125
  • Sep 18 06:59

    apbraun on exp-fix

    Changed the version number (compare)

  • Sep 17 20:37
    Datseris commented #125
  • Sep 17 18:15
    apbraun commented #125
  • Sep 17 15:28
    Datseris commented #125
  • Sep 17 15:01
    apbraun opened #125
  • Sep 17 14:55

    apbraun on exp-fix

    Changed `e` to `Base.MathConsta… (compare)

  • Sep 17 14:55

    apbraun on exp-fix

    (compare)

  • Sep 12 17:22

    Datseris on gh-pages

    build based on 1c386a2 (compare)

  • Sep 12 17:11

    github-actions[bot] on v1.5.0

    (compare)

  • Sep 12 16:08

    Datseris on gh-pages

    build based on 1c386a2 (compare)

  • Sep 12 15:57
    JuliaRegistrator commented on 1c386a2
Farhan Omar
@FO2020
The code for the trace
%matplotlib inline  
#%matplotlib notebook
#%pylab
from matplotlib import use
#use("Qt5Agg")
import numpy as np
from numpy import *
from cmath import*
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.signal import find_peaks
from scipy import signal
import scipy as scp
from scipy.integrate import solve_ivp, odeint
import pandas as pd 
from scipy import linalg as LA
import nolds

chaotic=np.array([3.99256592e-04, 2.80065918e-01, 9.78851318e+00, 4.33094482e+00,
       1.59216309e-01, 4.10463217e+00, 7.37010376e-01, 1.67592651e-01,
       4.92309570e+03, 7.03787720e+02, 9.29061279e+01])
periodic=np.array([2.02604065e-04, 1.64956665e-01, 8.81476746e+00, 3.22366028e+00,
        3.05567932e+00, 4.00346842e+00, 1.00645364e+00, 8.62795715e-02,
        6.99554443e+03, 8.34424774e+02, 4.04785461e+01])
# Model ,functions and solvers 

def DM_bacteria(t,X,pm):    
    x, y1, y2 = X
    a, b,delta1,delta2,psi,phi,omega1,omega2=pm
    dxdt=x*(1-a*x-y1-b*y2)
    dy1dt=y1*(-delta1-omega1*y1-psi*b*y2+x)
    dy2dt=y2*(-delta2-omega2*b*y2+phi*y1+b*x)
    return [dxdt,dy1dt,dy2dt]
def DM_jac(t,X,pm):
    x, y1,y2 = X
    a, b,delta1,delta2,psi,phi,omega1,omega2=pm
    J=np.array([[-2*a*x - b*y2 - y1 + 1, -x, -b*x], 
                [y1, -b*psi*y2 - delta1 - 2*omega1*y1 + x, -b*psi*y1], 
                [b*y2, phi*y2, -2*b*omega2*y2 + b*x - delta2 + phi*y1]])
    return J

# Calculate the trace of the Jacobian matrix for an orbit
#to check if the system is dissapative (  If  the trace is negative, then the system is dissapative)

def MeanTraceJac(param_values,IC,tMax,solvers,RelTol,AbsTol):
    a, b,delta1,delta2,psi,phi,omega1,omega2,x0,y10,y20,*_ = param_values
    t=[0,tMax]
    t_eval_pts=np.linspace(0,tMax,300*tMax)

    pm=np.array([ a,b,delta1,delta2,psi,phi,omega1,omega2])
    init=IC
    sol=solve_ivp(DM_bacteria,t,init,method= solvers,t_eval= t_eval_pts,args=(pm,),jac=DM_jac,rtol=RelTol,atol=AbsTol);
    x,y1,y2=sol.y
    TrJAll=[]
    for i,j,k in zip(x,y1,y2):
        x=i
        y1=j
        y2=k
        TrJ=-2*a*x - 2*b*omega2*y2 - b*psi*y2 + b*x - b*y2 - delta1 - delta2 - 2*omega1*y1 + phi*y1 + x - y1 + 1
        TrJAll.append(TrJ)
    TrJMean=np.mean(np.array(TrJAll))
    return TrJMean, TrJAll[-100:],sol.y


JJ1=MeanTraceJac(chaotic,np.array([10.2 , 0.416, 2.207]),10000,'RK45',1e-8,1e-8)
JJ2=MeanTraceJac(periodic,np.array([11.155,  0.391,  3.893]),10000,'RK45',1e-8,1e-8)
Farhan Omar
@FO2020
Julia code used in the calculation of LEs
#Import packages
using DifferentialEquations
using DynamicalSystems
using Plots
chaotic=Array([3.99256592e-04, 2.80065918e-01, 9.78851318e+00, 4.33094482e+00,
       1.59216309e-01, 4.10463217e+00, 7.37010376e-01, 1.67592651e-01])
periodic4=[2.02604065e-04, 1.64956665e-01, 8.81476746e+00, 3.22366028e+00,
        3.05567932e+00, 4.00346842e+00, 1.00645364e+00, 8.62795715e-02]
# Model ,functions and solvers

function DM_bacteria(dXdt,X,pm,t)    
    x, y1, y2 = X
    dxdt,dy1dt,dy2dt=dXdt
    a, b,delta1,delta2,psi,phi,omega1,omega2=pm
    dXdt[1]=x*(1-a*x-y1-b*y2)
     dXdt[2]=y1*(-delta1-omega1*y1-psi*b*y2+x)
     dXdt[3]=y2*(-delta2-omega2*b*y2+phi*y1+b*x)

end

function DM_jac(J,X,pm,t)
    x, y1,y2 = X
    a, b,delta1,delta2,psi,phi,omega1,omega2=pm
    J[1,1]=-2*a*x - b*y2 - y1 + 1
    J[2,1]=y1
    J[3,1]=b*y2

     J[1,2]= -x
     J[2,2]=  -b*psi*y2 - delta1 - 2*omega1*y1 + x
     J[3,2]= phi*y2

     J[1,3]= -b*x
     J[2,3]= -b*psi*y1
     J[3,3]=-2*b*omega2*y2 + b*x - delta2 + phi*y1

end

X01=[10.2  ;  0.416;  2.207]
tspan1 = (0.0,1000.0)
f = ODEFunction(DM_bacteria, jac=DM_jac)
pb_jac1 = ODEProblem(f,X0,tspan1,chaotic)

tspan2 = (0.0,1000.0)
X02=[11.155,  0.391,  3.893]
pb_jac2 = ODEProblem(f,X02,tspan2,periodic4)
SOL2 = solve(pb_jac2,Tsit5(),reltol=1e-8,abstol=1e-8,saveat=0.001)

ds1 = ContinuousDynamicalSystem(DM_bacteria, X01, chaotic, DM_jac)
LEs1 = lyapunovs(ds1, 100000,Ttr = 2000)
ds2 = ContinuousDynamicalSystem(DM_bacteria, X02, periodic4, DM_jac)
LEs2 = lyapunovs(ds2, 100000,Ttr = 2000)
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> did you try running it without the jacobian just to double check that it wasn't a mis-specification in that portion?
Farhan Omar
@FO2020
I will try running it without Jacobian and see the result will changes
Farhan Omar
@FO2020
@GitterIRCbot I removed the jacobian from ds and ran the code but got the same results.
BridgingBot
@GitterIRCbot
[slack] <chrisrackauckas> okay, and you've confirmed the ODE is correct?
Farhan Omar
@FO2020
@GitterIRCbot I run the solution via python , Matlab and plotted the solution and it is periodic . I have also analytical expression for the equilibrium points and all of them are unstable so the largest lyapunov exponent can't be negative for parameter set periodic 4
BridgingBot
@GitterIRCbot
[slack] <SebastianM-C> You could try using the lyapunov convergence function to investigate what happens during the calculation.
BridgingBot
@GitterIRCbot

[slack] <SebastianM-C> You could also see if the solver and it's parameters influence the results. I would try with Vern9. For the lyapunovs function, I think you can pass the solver via a keyword argument, but I don't recall it exactly.

And you can also try the largest lyapunov exponent function, which uses a different algorithm

Farhan Omar
@FO2020
I tried lyapunov also and the results matched with lyapunovs
@GitterIRCbot I will check with your suggestion and see. thanks
BridgingBot
@GitterIRCbot
[slack] <Datseris> Yeap, thanks for pointing this out... @asinghvi17 I think you contributed the eachcol method extention. We should have changed the limits in that PR...
[slack] <Datseris> Hi,
[slack] <Datseris> @FO2020, see this page for making a DynamicalSystem: https://juliadynamics.github.io/DynamicalSystems.jl/dev/ds/general/
[slack] <Datseris> Sebastian's suggestions are all great, and I've got nothing more to add: use higher-accuracy solvers, higher tolerances, as well as more integration time.
BridgingBot
@GitterIRCbot
[slack] <Datseris> You can do this with a code like
using OrdinaryDiffEq # access advanced solvers lyapunovs(ds, 100_000, k; alg = Vern9(), reltol = 1e-9, abstol = 1e-9)
(or any other args)
Farhan Omar
@FO2020
The results has improved after increasing the time to ten million and changed the solver as you suggested. Here are the results for periodic4
ds2 = ContinuousDynamicalSystem(DM_bacteria, X02, periodic4, DM_jac)
LEs2=lyapunovs(ds2, 10000000,Ttr = 5000, lg = Vern9(), reltol = 1e-10, abstol = 1e-10)
3-element Array{Float64,1}:
 -0.0001959745405066824
 -0.20996471314933712
 -0.20996467013046127
BridgingBot
@GitterIRCbot
[slack] <Datseris> great, thanks a lot to @SebastianM-C!
[slack] <Datseris> Also FO2020 the function lyapunov uses a different algorithm to get the maximum Lyapunov exponent (which in your case is enough to show whether you have distinguished between periodic and chaotic trajectory). You can always use both lyapunov and lyapunovs to be more certain of your results.
Farhan Omar
@FO2020
@GitterIRCbot Yes it is enough with LLE but I wanted all to check them with the sum of the trace . Thanks for the kind help
Mesut Karakoç
@mkarakoc
Hello to the community. I have just registered today! I wonder if it is possible to create obstacles in DynamicalBillards.jl which their shape/size can change with time? My aim is to simulate bouncing balls on a vibrating plate. If DynamicBillards.jl is not suitable for this purpose. Is there any other Library which it can achieve this?
BridgingBot
@GitterIRCbot
[slack] <Datseris> Hi mkarakoc, I am very interested in what you are saying. I also want to make DynamicalBilliards.jl support obstacles whose position and/or shape changes with time
[slack] <Datseris> About 2 months ago I wanted to make an animation of the stadium slowly changing into the circle billiard, to highlight the fundamental difference between chaotic and quasi-periodic trajectories.
[slack] <Datseris> (of course, the functionality is not there yet)
[slack] <Datseris> Here is a suggestion: why don't you open an issue at DynamicalBilliards.jl on github, and we discuss there more details, to find out the simpler way for this to happen?
[slack] <Datseris> I do have to say though that particles in DynamicalBilliards.jl are point particles, and if you want to have real balls bouncing around (with size), then you would need a new particle type (for straight propagation, I don't think this is too hard to do)
BridgingBot
@GitterIRCbot
[slack] <dpsanders> For real balls you can use HardSphereDynamics
[slack] <dpsanders> But then it’s hard to have different table shapes. So somehow we need to combine the two pacakges
[slack] <Datseris> I totally forgot to mention HardSphereDynamics! I always assosiate it as being able to do particle interaction (muuuuch more advanced and difficult concept), so I forgot that this was an option (since, you know, you can evolve a single particle)
BridgingBot
@GitterIRCbot
[slack] <dpsanders> You can also have “infinite-mass” particles that are fixed. You literally set their mass to Inf
[slack] <dpsanders> There’s a clever trick in the implementation of the collision rules to allow that to work 😉
[slack] <dpsanders> (syntax may be out of date?)
Mesut Karakoç
@mkarakoc
I will look into detail. Thank you so much for the help. I am a Pythonian but considering to use Julia. The ball in consideration is a point particle for now.
The problem I have interested is given in the following: https://journals.aps.org/pre/abstract/10.1103/PhysRevE.79.055201
BridgingBot
@GitterIRCbot
[slack] <dpsanders> I’m not sure how useful the package is for that use case. Do they describe their algorithm? You should just be able to implement it directly using some kind of root finder
[slack] <dpsanders> Yes they basically do. You can use e.g. Roots.jl or IntervalRootFinding.jl
[slack] <dpsanders> (weird and interesting paper)
BridgingBot
@GitterIRCbot
[slack] <Datseris> I agree with David for this specific case. Regardless, having the possibility of obstacles being able to move would be great for DynamicalBilliards.jl, so I'd be happy to guide someone that wants to do such a PR. It would also be great to combine functionality, or bridge usage between HardSphereDynamics and DynamicalBilliards, but I feel this is not as simple as a PR.
Marcelo Forets
@mforets
Hi @Datseris ! To expand a bit on my discourse answer (here),
Screenshot from 2020-09-03 23-57-59.png
the time interval associated to the sets in red gives both a lower bound on the time entering the ball, and an upper bound on the time exiting the ball.
Mesut Karakoç
@mkarakoc

@GitterIRCbot I have obtained their results using my own algorithm (written as a Python code) with the Newton-Raphson method. I had/have a feeling DynamicalBilliards.jl (or some other package of JuliaDynamics) kind of making simulations (!?) rather than calculations of the trajectory equations. I might be wrong about my feelings :) . So, if i can make some simulations (kind of experiments in a simulated world), i might be able to compare my own results with the simulations (experiments).

The CAD software in the video might be another option: https://www.youtube.com/watch?v=tMKXbLBgkEc&t=635s
I also thought some game engines but they will be over-complicated to learn for now.

BridgingBot
@GitterIRCbot
[slack] <Datseris> Haha, that is a great video, but not the kind of simulations targeted by JuliaDynamics, at least to my knowledge. Perhaps something in JuliaRobotics might be more fitting?
BridgingBot
@GitterIRCbot
[slack] <dpsanders> I don’t understand the distinction that you’re making between “simulations” and “calculations of the trajectory equations”. To me they seem like the same thing
[slack] <Datseris> Yeah to me too, but perhaps this is some kind of different lingo from different schools of thought... so I thought of not asking 😛
[slack] <Datseris> The cited video deals more with material simulations than dynamics
Mesut Karakoç
@mkarakoc
I don't understand my self about the distinction, neither. :) Could a scientist (theoretician) in "The Matrix" (https://en.wikipedia.org/wiki/The_Matrix) make calculations and compare with the experimental results which are obtained by an experimentalist in The Matrix?
BridgingBot
@GitterIRCbot
[slack] <jamesjscully> These simulations are essentially boundary value problems for PDEs so they are dynamical systems. You can theoretically set up the problem in DifferentialEquations.jl, but with CAD software it's MUCH easier to configure