Cos'è ElectroYou | Login Iscriviti

ElectroYou - la comunità dei professionisti del mondo elettrico

Ricerca personalizzata
11
voti

Simbolica o Numerica ?

Indice

Abstract

Soluzione di un'equazione differenziale nella modalita' simbolica e numerica.


Introduzione

Prendendo spunto da un post del forum di elettronica su un circuito R-L-C in regime transitorio, analizziamo la possibilita' di risoluzione al calcolatore di una equazione differenziale del secondo ordine, attraverso l'uso di tre programmi gratuiti presenti nell'elenco dei Free Tools di ElectroYou.

Per la soluzione simbolica useremo Maxima, mentre per la soluzione numerica useremo VisSim e Scilab.

Circuito

Dato il seguente circuito R-L-C

fig.0

fig.0

alimentato da un generatore di corrente definito dalle seguenti relazioni


\left\{ \begin{align}
  & J(t)=-1\,A\,\,\,\,t<0 \\ 
 & J(t)=+1\,A\,\,\,\,t\ge 0 \\ 
\end{align} \right.


dal principio di Kirchhoff al nodo A otteniamo


iR + iL + iC = 1


e applicando le equazioni caratteristiche di ogni bipolo avremo


\frac{v}{R}+\frac{1}{L}\int{vdt+}\,C\frac{dv}{dt}=1


che per comodita' riscriviamo come


\frac{1}{RC}\cdot \frac{dv}{dt}+\frac{1}{LC}\cdot v\,+\frac{d^{2}v}{dt^{2}}=0


Ci proponiamo di ricercare una soluzione per via simbolica con wxMaxima e una soluzione per via numerica con VisSim e Scilab

Maxima

Per la simulazione per via simbolica, facciamo uso di xwMaxima (purtroppo Maple Demo non permette la risoluzione di equazioni differenziali).

Inseriamo l'equazione nella finestra principale come prima "cella" ricordando di terminare con un ";" e completare con un SHIFT+ENTER .

fig.1

fig.1

Notiamo come automaticamente Maxima nomini la riga di input con (%i30) e la riga di output con (%o30).

La successiva richiesta di risoluzione della ODE (Ordinary Differential Equation) di secondo ordine (2) avviene alla riga (%i31) con un ode2; il "%" indica l'ultimo risultato dell'elaborazione (equivalente all' "ans" di Scilab).

Maxima a questo punto, non avendo ancora assegnato i valori numerici di R, L e C, richiede per proseguire, l'input del segno di L\left( 4CR^{2}-L \right) che corrisponde curiosamente, non al segno del discriminante Δ, ma del suo opposto; Maxima prende cioe' come scelta principale il caso di radici complesse, nel nostro esempio avendo supposto per semplicita' R=C=L=1, la nostra risposta sara' "p" (lascio al lettore la simulazione con le altre due alternative).

Ricordo che il discriminante di un'equazione di secondo grado e' pari a


\begin{align}
  & \Delta =b^{2}-4ac \\ 
 & \Delta =\frac{1}{R^{2}C^{2}}-\frac{4}{CL}=\frac{L-4R^{2}C}{R^{2}C^{2}L}=\frac{L\left( L-4R^{2}C \right)}{R^{2}C^{2}L^{2}} \\ 
\end{align}


L'output (%o31) e' la soluzione generale dell'equazione differenziale; notiamo la presenza delle costanti %K1 e %K2, dei due termini in seno e coseno (caratteristici delle radici complesse) e del termine di smorzamento e − βt.

Con (%i32) inseriamo la riga delle condizioni iniziali attraverso ic (initial conditions) considerando i valori iniziali al tempo t=0+.

Non essendo ammesse discontinuita' sulla tensione ai capi del condensatore e sulla corrente nell'induttore,


\frac{v(0+)}{R}+i_{L}(0+)+C\frac{dv}{dt}(0+)=1\,\,\,\,\to \,\,\frac{0}{R}-1+C\frac{dv}{dt}(0+)=1


\frac{dv}{dt}(0+)=\frac{2}{C}


per mantenere una maggior generalita' inseriamo per t=0+, v=v0 e dv/dt=2/C e otteniamo l'equazione finale (%032) per v(t).

E' interessante notare come, con un semplice comando tex(%), possiamo ottenere la conversione della formula in codice TEX;

 tex(%);
 $$v=\left({{\sin \left({{t\,\sqrt{{{4}\over{C\,L}}-{{C^2}\over{R^2}}}
  }\over{2}}\right)\,\left(4\,R+{\it v_0}\,C^2\right)}\over{C\,\sqrt{
 {{4}\over{C\,L}}-{{C^2}\over{R^2}}}\,R}}+{\it v_0}\,\cos \left({{t\,
 \sqrt{{{4}\over{C\,L}}-{{C^2}\over{R^2}}}}\over{2}}\right)\right)\,e
 ^ {- {{t\,C}\over{2\,R}} }$$

un semplice copia-incolla, ci permetterà di scriverla, in formato HQ;


(%032)  v=\left({{\sin \left({{t\,\sqrt{{{4}\over{C\,L}}-{{C^2}\over{R^2}}}
 }\over{2}}\right)\,\left(4\,R+{\it v_0}\,C^2\right)}\over{C\,\sqrt{
 {{4}\over{C\,L}}-{{C^2}\over{R^2}}}\,R}}+{\it v_0}\,\cos \left({{t\,
 \sqrt{{{4}\over{C\,L}}-{{C^2}\over{R^2}}}}\over{2}}\right)\right)\,e
 ^ {- {{t\,C}\over{2\,R}} } .


A questo punto (%i34), particolarizziamo v(t) fissando le costanti del circuito: R=1, L=1, C=1, V0=0 ottenendo la (%o34), che rappresenta sia la tensione v(t) sul parallelo sia la corrente IR(t).


fig.2

fig.2


Per il calcolo della corrente nel condensatore IC(t) e di quella nell'induttore IL(t), passiamo a differenziare (%i46) ed a integrare (%i88) la tensione sul condensatore v(t).


Con l'integrazione, IL(t) e' determinata a meno di una costante e di conseguenza valutiamo x per t=0 (%i92); il risultato pari a -2 ci permette di determinare la costante.

Da IL(0)+k=-2+k con IL(0)=-1, ricaviamo k=1.


fig.3

fig.3

Abbiamo cosi' i risultati finali per le tre correnti IR(t)>(%o98), IC(t)>(%o99) e IL(t)>(%o100).


Per la visualizzazione grafica finale usiamo un plot2d e Maxima attraverso Gnuplot ci dara' il seguente grafico; si nota, a conferma della soluzione, come la somma delle tre correnti sia sempre pari a 1 ampere.

fig.4

fig.4

VisSim

La simulazione numerica attraverso VisSim si ottiene come al solito costruendo per prima cosa la catena di integrazione a partire dalla derivata prima dv/dt > funzione v(t) > integrale di v(t) [§v].


La simulazione del generatore di corrente con intensita' che passa da -1A a +1A al tempo zero, e' ottenuta attraverso un blocco sommatore che addiziona una costante pari a -1 ad uno scalino di ampiezza pari a 2 al tempo zero e le condizioni iniziali degli integratori entrambe a zero.

fig.5

fig.5

la struttura dello schema a blocchi e' costruita sulla base della seguente rappresentazione analitica


\frac{dv}{dt}=\frac{1}{C}\left[ J-\left( \frac{1}{L}\int{vdt}+\frac{v}{R} \right) \right]


dove l'uguaglianza fre i due membri dell'equazione differenziale e' rappresentata dall'ultimo collegamento di controreazione (in alto a sinistra).

Il collegamento dei segnali principali al display, e i settaggi sotto Simulate > Simulation Properties > Start=-10s, Time-Step=0.01s, End=10s, ed un GO finale, completano la simulazione.


VisSim permette di effettuare una simulazione dinamica; sostituendo le costanti di ingresso con degli "slider" viene permessa una regolazione continua dei parametri in un predeterminato campo di valori.

fig.5b

fig.5b

In questa seconda simulazione, a differenza della prima (fig.5), lo scalino di ingresso a J e' di 1 ampere mentre la condizione iniziale del secondo integratore e' stata impostata a -L.

[Errata Corrige: nei grafici di fig.5 e fig.5b sostituire (sec) con (s)]

Scilab

In Scilab il discorso e' simile a quello visto per VisSim anche se un po' piu' complesso dal punto di vista matematico.

Usiamo l'algoritmo associato a ode che rappresenta la funzione standard per la risoluzione delle equazioni differenziali del tipo


\frac{dy}{dt}=f(t,y)\,\,\,\,,\,\,\,\,y(0)=y_0


Per i gradi superiori al primo si e' costretti ad un'impostazione vettoriale nella quale le successive derivate vengono ricavate, per via numerica, con una serie di equazioni differenziali del primo ordine.


Il codice relativo sara':

   // **************************************
  // Equazione differenziale secondo ordine
 // y + y'/(RC) + y/(LC) = 0
// **************************************
   L=1;C=1;R=1;
// Definisco la funzione
   function [z]=diffeq(t,y)
     z(1)=y(2);
     z(2)=-y(2)/(R*C)-y(1)/(L*C) ;
   endfunction
// Condizioni iniziali
     y(2)=2/C; y(1)=0; 
     y0=[y(1);y(2)];
// Definiamo Start - Step - End 
     t0=0; t=0:.02:10; // 500 passi  
// Soluzione primo ordine
     y=ode(y0,t0,t,diffeq);
// Grafico della iR(t)
    subplot(3,1,1)
    plot2d(t,y(1,:),5);xgrid()
    xtitle("Correnti","t","i R (t)");
// Grafico della iC(t)
    subplot(3,1,2)
    plot2d(t,y(2,:),2);xgrid()
    xtitle(" ","t" ,"i C (t)");
// Grafico della iL(t)
   subplot(3,1,3)
   plot2d(t,1-y(1,:)-y(2,:),3);xgrid()
   xtitle(" ","t" ,"i L (t)");
//------------------------------------eof


Con l'output grafico per le correnti iR(t), iC(t) e iL(t) di fig.6.


fig.6

fig.6


Ed infine, costruiamo in modo simile a quanto fatto per VisSim, uno schema a blocchi con Scicos

che fornirà i seguenti andamenti per IL, IC e IR.

Ovviamente la costruzione dello schema a blocchi con Scicos può essere, come in questo caso, complessa; ci si può allora chiedere se esista un metodo più rapido per inserire il sistema lineare in Scicos.

Ricordando che in generale sarà possibile rappresentare il sistema attraverso le seguenti relazioni


\left\{ \begin{align}
  & \frac{dx_{1}}{dt}=a_{11}x_{1}+a_{12}x_{2}+b_{1}u \\ 
 & \frac{dx_{2}}{dt}=a_{21}x_{1}+a_{22}x_{2}+b_{2}u \\ 
 & y=c_{1}x_{1}+c_{2}x_{2}+du \\ 
\end{align} \right.


e dalle condizioni iniziali


\left\{ \begin{align}
  & x_{1}(0)=x_{10} \\ 
 & x_{2}(0)=x_{20} \\ 
\end{align} \right.


particolarizzando per il nostro sistema dinamico


\left\{ \begin{align}
  & \frac{dv_{c}}{dt}=\frac{1}{C}(J-i_{L}-i_{R}) \\ 
 & \frac{di_{L}}{dt}=\frac{v_{c}}{L} \\ 
 & i_{R}=\frac{v_{c}}{R} \\ 
\end{align} \right.


e ponendo


u=J,\,\,\,x_{1}=v_{c},\,\,\,x_{2}=i_{L},\,\,\,y=i_{R}


potremo scrivere il sistema di equazioni in forma matriciale come


\left\{ \begin{align}
  & \frac{dx}{dt}=Ax+Bu \\ 
 & y=Cx+Du \\ 
\end{align} \right.


con


A=\left[ \begin{matrix}
   -1/RC & -1/C  \\
   1/L & 0  \\
\end{matrix} \right]\,\,\,\,\,\,\,\,B=\left[ \begin{matrix}
   1/C  \\
   0  \\
\end{matrix} \right]\,\,\,\,\,\,\,C=\left[ \begin{matrix}
   1/R & 0  \\
\end{matrix} \right]\,\,\,\,\,\,D=0

e condizioni iniziali x0 = [0, − 1], tensione nulla e corrente pari a -1 ampere nell'induttore, usando l'Editor inseriremo

e quindi in Scicos potremo fare riferimento alle matrici A,B,C,D e x0 per costruire il blocco rappresentativo il sistema

settando i parametri con

e ottenendo dalla simulazione l'andamento della iR(t)

NB si fa notare come aggiungendo altri 2 blocchi CLSS in parallelo, si sarebbero ottenuti anche gli andamenti delle altre due correnti !


to be continued ...



1

Commenti e note

Inserisci un commento

di ,

Beh, non c'è che dire RenzoDF: a chi ha voglia di imparare sai offrire parecchie e sempre nuove opportunità.

Rispondi

Inserisci un commento

Per inserire commenti è necessario iscriversi ad ElectroYou. Se sei già iscritto, effettua il login.