# Simulation of linear circuits

## 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
```

```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
```

```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