Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

About chua circuit numerical solutions

  1. Dec 1, 2016 #1
    hi all !!
    i'm trying to solve numerically the chua circuit set of differential equations , im using the equations showed in this pdf:
    http://nonlinear.eecs.berkeley.edu/chaos/RobustOpAmpRealizationOfChuaCircuit.pdf

    i have the real circuit mounted and i'm using its parameters for the numerical solution

    i'm using the RK4 method solving this system , but i got a really big problem calculating the k slopes :

    these all grow really fast without control!! i found that the problem is caused by the term (RC_i)^{-1} i=1,2 since the capacitances used are really small compared to the resistance used

    this is my awk code for solving this problem:

    awk '
    function abs(value){
    return (value<0?-value:value)
    }
    function g(MN1,M0,M1,BPK1,BKP2,VOLT1){
    return MN1*VOLT1+0.5*(M0-MN1)*(abs(VOLT1+BKP1)-abs(VOLT1-BKP1))+0.5*(M1-M0)*(abs(VOLT1+BKP2)-abs(VOLT1-BKP2))
    }
    function tvc1(RE,C1,MN1,M0,M1,BPK1,BKP2,VOLT1,VOLT2){
    return (1/(RE*C1))*(VOLT2-VOLT1-RE*g(MN1,M0,M1,BPK1,BKP2,VOLT1))
    }
    function tvc2(RE,C2,VOLT1,VOLT2,INTL){
    return (1/(RE*C2))*(VOLT1-VOLT2+RE*INTL)
    }
    function til(LE,VOLT2){
    return -VOLT2/LE
    }
    {}
    END{
    #variable parameters of the circuit
    R=1.2e3
    RL=2.14e3
    #fixed parameters in the circuit
    r1=220.
    r3=2.2e3
    r4=22.0e3
    r6=3.3e3
    c1=10.e-9
    c2=100.e-9
    #equivalent inductance given by gyrator circuit
    r8=1.e3
    r9=1.e3
    r10=100.
    c3=100.e-9

    L=(r10*r8*c3*RL)/r9

    #parameters for the function of chua diode
    Esat=8.3
    mn1=1/r1+1/r4
    m0=1/r4-1/r3
    m1=-1/r6-1/r3
    bp1=r3*Esat/(r3+r1)
    bp2=r6*Esat/(r6+r4)

    print mn1,m0,m1,bp1,bp2
    #USING THE RK4 METHOD FOR SOLVING THE SYSTEM OF DIFFERENTIAL EQUATIONS
    #real life time of circuit simulation
    T=0.01
    #number of temporal points
    n=1000.
    #size of temporal points
    h=T/n
    #initial time
    t=0.
    #initial conditions
    v1=3.
    v2=1.
    il=1.
    print v1,v2,il
    #inicializating variables
    vn1=0.
    vn2=0.
    iln=0.
    print (1/(R*c1))*(v2-v1-R*g(mn1,m0,m1,bp1,bp2,v1))
    while(t<T){
    #finding slopes with the k's
    # k1
    kvn11=tvc1(R,c1,mn1,m0,m1,bp1,bp2,v1,v2)
    kvn21=tvc2(R,c2,v1,v2,il)
    kiln1=til(L,v2)


    # k2
    kvn12=tvc1(R,c1,mn1,m0,m1,bp1,bp2,v1+kvn11*h/2.,v2+kvn21*h/2.)
    kvn22=tvc2(R,c2,v1+kvn11*h/2.,v2+kvn21*h/2.,il+kiln1*h/2.)
    kiln2=til(L,v2+kvn21*h/2.)
    # k3
    kvn13=tvc1(R,c1,mn1,m0,m1,bp1,bp2,v1+kvn12*h/2.,v2+kvn22*h/2.)
    kvn23=tvc2(R,c2,v1+kvn12*h/2.,v2+kvn22*h/2.,il+kiln2*h/2.)
    kiln3=til(L,v2+kvn22*h/2.)
    # k4
    kvn14=tvc1(R,c1,mn1,m0,m1,bp1,bp2,v1+kvn13*h,v2+kvn23*h)
    kvn24=tvc2(R,c2,v1+kvn13*h,v2+kvn23*h,il+kiln3*h)
    kiln4=til(L,v2+kvn23*h)

    #finding new points
    vn1=v1+(1./6.)*(kvn11+2*kvn12+2*kvn13+kvn14)*h
    vn2=v2+(1./6.)*(kvn21+2*kvn22+2*kvn23+kvn24)*h
    iln=il+(1./6.)*(kiln1+2*kiln2+2*kiln3+kiln4)*h
    #terminal print of new set of points
    print vn1,vn2,iln
    #next iteration preparation
    v1=vn1
    v2=vn2
    il=iln
    t=t+h
    }
    }
    ' $1

    any suggestion for the code ? or to solve my problem of excessive slope grow ?
    thanks for reading my post :)
     
  2. jcsd
  3. Dec 7, 2016 #2
    Thanks for the thread! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post? The more details the better.
     
  4. Dec 9, 2016 #3

    Baluncore

    User Avatar
    Science Advisor

    BKP1 or BPK1 in function g() definition ?
    function g(MN1,M0,M1,BPK1,BKP2,VOLT1){
    function g(MN1,M0,M1,BKP1,BKP2,VOLT1){
     
  5. Dec 10, 2016 #4

    Baluncore

    User Avatar
    Science Advisor

    Are awk variable names case sensitive ?

    I think h needs to be closer to 1e-6
    Start with T = 0.001 and t = 1000
    You must tune the circuit. Change R to about 2000.
    Use x and y range from -10 to +10. The plot of V1 against V2 shows the double scroll.
     
    Last edited: Dec 10, 2016
  6. Feb 22, 2017 #5

    Baluncore

    User Avatar
    Science Advisor

    Attached is the LTspice model and plot file. To run it, remove the .txt extensions.
    Likewise rename the file extension to run the Freebasic code.
     

    Attached Files:

Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: About chua circuit numerical solutions
  1. Numerical methods (Replies: 1)

  2. Numerical analysis (Replies: 1)

Loading...