- #1

chastiell

- 11

- 0

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 :)