About chua circuit numerical solutions

[chastiell]

hi all !
i'm trying to solve numerically the chua circuit set of differential equations , I am 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 :)

 

[Baluncore]

BKP1 or BPK1 in function g() definition ?
function g(MN1,M0,M1,BPK1,BKP2,VOLT1){
function g(MN1,M0,M1,BKP1,BKP2,VOLT1){

 

[Baluncore]

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.

 

[Baluncore]

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.