1. Resolution method
1.a. Laplace transmittance
We consider a linear circuit comprising N nodes and we denote by Vi the electric potential of node i. The circuit comprises an ideal voltage source (the input of the filter) whose potential will be taken by convention equal to 1. By definition, the ground is a node of the circuit whose potential is zero.
For a node k to which passive linear dipoles are attached, the law of nodes is written: ∑iYiVk-∑iYiVi = 0
where the sum relates to the nodes connected to node k by an admittance dipole Yi. The admittances are defined as functions of the variable s = jω.
The system of linear equations for the potentials takes the form AV + B = 0 where A is an N × N matrix of rational fractions, V and B two column vectors. Row k of matrix A corresponds to node k.
The addition of a source of direct voltage (constant) between the nodes na and nb introduces for these two nodes the equation: Vna-Vnb = U
A voltage-controlled voltage source (STCT) between nodes na and nb, controlled by the d.d.p. between nodes nx and ny, introduce the equation: Vna-Vnb = g (Vnx-Vny)
for the nodes na and nb.
If the output of the filter is the voltage at node ns, the Laplace transmittance of the filter is the potential of that node: h (s) = Vns (s)
1.b. Frequency response
The response in sinuoidal regime is given by the transfer function: H (ω) = h (jω)
1.c. Temporal response
It is assumed that at the initial time t = 0 the capacitors are discharged and the inductors have no current. In this case, the Laplace transform of the output signal is: L (s) = h (s) F (s)
where F (s) is the Laplace transform of the input signal. For a voltage step applied at the input, we have: F (s) = 1s
The impulse response is obtained with: F (s) = 1
2. Simulation with Scilab
2.a. System definition
The following functions transform either the matrix A only, or the matrixes A and B. The voltage sources must be defined after the passive dipoles because they deactivate the law of the nodes.
function [P]=ajouter_dipole(A,na,nb,Y)
P=A
P(na,na)=P(na,na)+Y
P(na,nb)=P(na,nb)-Y
P(nb,nb)=P(nb,nb)+Y
P(nb,na)=P(nb,na)-Y
endfunction
Addition of special passive dipoles:
function [P]=ajouter_resistance(A,na,nb,R)
P=ajouter_dipole(A,na,nb,1/R)
endfunction
function [P]=ajouter_capacite(A,na,nb,C)
P=ajouter_dipole(A,na,nb,poly(0,'s')*C)
endfunction
function [P]=ajouter_inductance(A,na,nb,L)
P=ajouter_dipole(A,na,nb,1/(poly(0,'s')*L))
endfunction
Adding a DC voltage source:
function [P,Q]=ajouter_source_tension_continue(A,B,na,nb,U)
P=A
Q=B
N=size(P,'c')
P(na,:)=zeros(1,N)
P(nb,:)=zeros(1,N)
P(na,na)=1
P(na,nb)=-1
Q(na)=-U
P(nb,nb)=-1
P(nb,na)=1
Q(nb)=-U
endfunction
Addition of a voltage controlled voltage source:
function [P,Q]=ajouter_source_STCT(A,B,na,nb,nca,ncb,g)
P=A
Q=B
N=size(P,'c')
P(na,:)=zeros(1,N)
P(nb,:)=zeros(1,N)
P(na,na)=1
P(na,nb)=-1
P(na,nca)=-g
P(na,ncb)=g
Q(na)=0
P(nb,nb)=-1
P(nb,na)=1
P(nb,nca)=-g
P(nb,ncb)=g
Q(nb)=0
endfunction
Add a mass, i.e. cancel a potential:
function [P,Q]=ajouter_masse(A,B,n0)
P=A
Q=B
P(n0,:)=zeros(1,size(P,'c'))
P(n0,n0)=1
Q(n0)=0
endfunction
If ne designates the node of the input of the filter, we set Vne = 1. In this way, the potential Vk directly gives the transfer function obtained by placing the output on the node k:
function [P,Q]=definir_entree(A,B,ne)
P=A
Q=B
N=size(P,'c')
P(ne,:)=zeros(1,N)
P(ne,ne)=1
Q(ne)=-1
endfunction
2.b. Frequency response
For each pulse ω, the matrix of rational fractions A is evaluated with the horner function, then the system is solved.
function [w,db,phi]=reponse(A,B,noeud,lgmin,lgmax,np)
w=logspace(lgmin,lgmax,np)
H=w
for i=1:np,
M=horner(A,%i*w(i))
[x0,ker]=linsolve(M,B) H(i)=x0(noeud) end;
[db,phi]=dbphi(H) endfunction
2.c. Example
Consider the following RLC circuit:
The following function defines the circuit:
getf('simulineaire.sci');
function [A,B]=circuitRLC(R,L,C)
N=4
A=zeros(N,N)
B=zeros(N,1)
A=ajouter_resistance(A,1,2,R)
A=ajouter_inductance(A,2,3,L)
A=ajouter_capacite(A,3,4,C)
[A,B]=ajouter_masse(A,B,4)
[A,B]=definir_entree(A,B,1)
endfunction
Definition of the circuit and plot of the Bode diagram for node 3
[A,B]=circuitRLC(10,1e-3,1e-6);
[w,db,phi]=reponse(A,B,3,2,6,200);
plot1=scf();
plot2d(w,db,style=5,logflag='ln');
xtitle('Circuit RLC','log f','GdB');
plot2=scf();
plot2d(w,phi,style=5,logflag='ln');
xtitle('Circuit RLC','log f','phi');
3. Simulation with Mathematica
Mathematica module:
3.a. System definition
ajouterDipole[A_,na_,nb_,Y_]:=Module[{P},
P=A;
P[[na,na]]=P[[na,na]]+Y;
P[[na,nb]]=P[[na,nb]]-Y;
P[[nb,nb]]=P[[nb,nb]]+Y;
P[[nb,na]]=P[[nb,na]]-Y;
Return[P];
]
ajouterResistance[A_,na_,nb_,R_]:=ajouterDipole[A,na,nb,1/R];
ajouterCapacite[A_,na_,nb_,C_]:=ajouterDipole[A,na,nb,C*s];
ajouterInductance[A_,na_,nb_,L_]:=ajouterDipole[A,na,nb,1/(L*s)];
ajouterSourceTensionContinue[A_,B_,na_,nb_,U_]:=Module[{P,Q,d,n},
P=A;
Q=B;
d = Dimensions[A];
n = d[[2]];
P[[na]]=Table[0,{n}];
P[[nb]]=Table[0,{n}];
P[[na,na]]=1;
P[[na,nb]]=-1;
P[[nb,nb]]=-1;
P[[nb,na]]=1;
Q[[nb]]=-U;
Return[{P,Q}];
]
ajouterSourceTensionSTCT[A_,B_,na_,nb_,nx_,ny_,g_]:=Module[{P,Q,d,n},
P=A;
Q=B;
d = Dimensions[A];
n = d[[2]];
P[[na]]=Table[0,{n}];
P[[nb]]=Table[0,{n}];
P[[na,na]]=1;
P[[na,nb]]=-1;
P[[nb,nb]]=-1;
P[[nb,na]]=1;
P[[na,nx]]=-g;
P[[na,ny]]=g;
P[[nb,nx]]=-g;
P[[nb,ny]]=g;
Q[[nb]]=0;
Return[{P,Q}];
]
ajouterMasse[A_,B_,n0_]:=Module[{P,Q,d,n},
P=A;
Q=B;
d = Dimensions[A];
n = d[[2]];
P[[n0]]=Table[0,{n}];
P[[n0,n0]]=1;
Q[[n0]]=0;
Return[{P,Q}];
]
definirEntree[A_,B_,ne_]:=Module[{P,Q,d,n},
P=A;
Q=B;
d = Dimensions[A];
n = d[[2]];
P[[ne]]=Table[0,{n}];
P[[ne,ne]]=1;
Q[[ne]]=-1;
Return[{P,Q}];
]
3.b. Transmittance
Laplace transmittance for an output on node ns:
transfert[A_,B_,ns_]:=Module[{v},
v=LinearSolve[A,-B];
Return[v[[ns]]];
]
3.c. Frequency response
bodeGain[h_,lgmin_,lgmax_,Gmin_,Gmax_]:=Module[{},
Return[Plot[20*Log[10,Abs[h/.{s->I*2*Pi*10^p}]],{p,lgmin,lgmax},PlotRange->{{lgmin,lgmax},{Gmin,Gmax}},AxesLabel->{"log f","GdB"}]];
]
bodePhase[h_,lgmin_,lgmax_]:=Module[{},
Return[Plot[180/Pi*Arg[h/.{s->I*2*Pi*10^p}],{p,lgmin,lgmax},AxesLabel->{"log f","phi"}]];
]
3.d. Response to a step
The following function takes the Laplace transmittance and the maximum time as argument, and returns the time response curve:
reponseEchelon[h_,tmax_,vmin_,vmax_]:=Module[{r},
r=InverseLaplaceTransform[h/s,s,t];
Return[Plot[r/.{t->temps},{temps,0,tmax},PlotRange->{{0,tmax},{vmin,vmax}},AxesLabel->{"t","V"}]];
]
The impulse response is defined in the same way:
reponseImpulsion[h_,tmax_,vmin_,vmax_]:=Module[{r},
r=InverseLaplaceTransform[h,s,t];
Return[Plot[r/.{t->temps},{temps,0,tmax},PlotRange->{{0,tmax},{vmin,vmax}},AxesLabel->{"t","V"}]];
]
3.e. Example
We take the example discussed above:
Get['simulineaire.m'];
circuitRLC[R_,L_,C_]:=Module[{A,B,n},
n=4;
A=Table[0,{n},{n}];
B=Table[0,{n}];
A=ajouterResistance[A,1,2,R];
A=ajouterInductance[A,2,3,L];
A=ajouterCapacite[A,3,4,C];
{A,B}=ajouterMasse[A,B,4];
{A,B}=definirEntree[A,B,1];
Return[{A,B}];
]
{A,B}=circuitRLC[10,10^-3,10^-6];
Transfer function for node 3:
h=transfert[A,B,3]
1000000000/(1000000000 + 10000*s + s^2)
Diagramme de Bode :
bodeGain[h,2,6,-80,20]
bodePhase[h,2,6]
Response to a step:
reponseEchelon[h,2*10^-3,0,2]
The calculation can be done formally:
{A,B}=circuitRLC[R,L,C];
h=transfert[A,B,3]/.s->I*omega
1-C L ω2+i C ω R+1
