Diferencia entre revisiones de «Circuitos Eléctricos RL»

De MateWiki
Saltar a: navegación, buscar
(Método de Euler)
(Modificación de las intensidades)
 
(No se muestran 83 ediciones intermedias del mismo usuario)
Línea 15: Línea 15:
 
donde <math>L</math> es el coeficiente de autoinducción en henrios (<math>H</math>).
 
donde <math>L</math> es el coeficiente de autoinducción en henrios (<math>H</math>).
  
También tenemos en cuenta las '''Leyes de Kirchoff''', que establecen el comportamiento de los circuitos:
+
También tenemos en cuenta las '''Leyes de Kirchhoff''', que establecen el comportamiento de los circuitos:
 
# '''Ley de corrientes''': En cada nodo, la suma de corrientes que entra es igual a la que sale.
 
# '''Ley de corrientes''': En cada nodo, la suma de corrientes que entra es igual a la que sale.
 
# '''Ley de tensiones''': En cada ciclo cerrado o malla, la suma de diferencias de potencial es nula.
 
# '''Ley de tensiones''': En cada ciclo cerrado o malla, la suma de diferencias de potencial es nula.
Línea 22: Línea 22:
 
[[Archivo:Ecs1.png|400px|miniatura|centro|Circuitos Eléctricos RL]]
 
[[Archivo:Ecs1.png|400px|miniatura|centro|Circuitos Eléctricos RL]]
  
== Ley de Kirchoff de Voltaje o Tensiones ==
+
== Ley de Kirchhoff de Voltaje o Tensiones ==
  
 
La ley de voltaje de Kirchhoff indica que la suma de voltajes alrededor de una trayectoria o circuito cerrado debe ser cero.
 
La ley de voltaje de Kirchhoff indica que la suma de voltajes alrededor de una trayectoria o circuito cerrado debe ser cero.
Línea 41: Línea 41:
 
[[Archivo:cc.jpg|450px|miniatura|derecha|Gráfico representación analítica intensidad-tiempo.]]
 
[[Archivo:cc.jpg|450px|miniatura|derecha|Gráfico representación analítica intensidad-tiempo.]]
  
Suponiendo que en el instante t=0 el circuito pasa de estar abierto a cerrado, obtenemos la intensidad para cada instante t>0 teniendo en cuenta que la fuente de alimentación tiene un voltaje constante de E=20V, L=0.2H y R=5Ω. Para ello, aplicamos el método de Euler mediante un problema de valor inicial (P.V.I.) y lo comparamos con los resultados que hemos obtenido analíticamente.
+
Suponiendo que en el instante t=0 el circuito pasa de estar abierto a cerrado, obtenemos la intensidad para cada instante t>0 teniendo en cuenta que la fuente de alimentación tiene un voltaje constante de E=20V, L=0.2H y R=5Ω. El cálculo analítico de la intensidad quedará resuelto de la siguiente manera, mediante un Problema de Valor Inicial (P.V.I.) o de Cauchy:
 
   
 
   
  
Línea 47: Línea 47:
  
 
=== Método de Euler===
 
=== Método de Euler===
 +
 +
Aplicamos el método de Euler basado en la aproximación del valor de la función a la recta tangente en cada punto conocido, mediante un problema de valor inicial (P.V.I.) y lo comparamos con los resultados que hemos obtenido analíticamente.
  
 
[[Archivo:Graf2.png|450px|miniatura|derecha|Gráfico intensidad respecto del tiempo cuando el circuito pasa de estar abierto a cerrado aplicando el método numérico.]]
 
[[Archivo:Graf2.png|450px|miniatura|derecha|Gráfico intensidad respecto del tiempo cuando el circuito pasa de estar abierto a cerrado aplicando el método numérico.]]
  
 
{{matlab|codigo=
 
{{matlab|codigo=
%Datos del problema:
+
%Datos del problema:  
L=0.2; %Autoinductancia de la bobina  
+
L=0.2;%Autoinductancia de la bobina  
R=5; %Resistencia
+
R=5;%Resistencia
 
tau =L/R; %Definición
 
tau =L/R; %Definición
 
E=20; %Valor de la fem
 
E=20; %Valor de la fem
t0=0; %Tiempo
+
t0=0;%Tiempo
 
tN=5*tau; %Tiempo de carga
 
tN=5*tau; %Tiempo de carga
 
i0=0; %En t=0 está descargado
 
i0=0; %En t=0 está descargado
Línea 63: Línea 65:
  
 
%Discretización:
 
%Discretización:
h=0.01;
+
h=0.001;
 
N=round((tN-t0)/h);
 
N=round((tN-t0)/h);
 
+
%Crear el vector t
%Crear el vector t:
+
 
t=linspace(t0,tN,N+1);
 
t=linspace(t0,tN,N+1);
 
i=zeros(1,N+1);
 
i=zeros(1,N+1);
 
i(1)=i0;
 
i(1)=i0;
 
+
%Euler
 
for k=1:N
 
for k=1:N
 
     i(k+1)=i(k)+h*f(t(k),i(k));
 
     i(k+1)=i(k)+h*f(t(k),i(k));
 
end
 
end
 +
ir=E/R*(1-exp(-t./tau)); %solción analítica
 +
clf
  
%Dibujamos:
+
%Representación:
plot(t,i)
+
hold on
 +
plot(t,i,'b')
 +
plot(t,ir,'r')
 +
xlabel('tiempo')
 +
ylabel('intensidad')
 +
legend('Euler','Exacta','Location','Best');
 +
hold off
 
}}
 
}}
  
Se ha elegido '''tN=5*tau=0.2s''' porque es el tiempo que tarda en cargarse la bobina.  
+
En la gráfica podemos ver que efectivamente el tiempo de carga es 5*tau.  
  
[[Archivo:TN3.jpg|450px|centro|Solución por el método de Euler para tN=3.]]
+
[[Archivo:H01_leyenda.png|550px|miniatura|centro|Solución por el método de Euler para h=0.01.]][[Archivo:H001_leyenda.png|550px|miniatura|centro|Solución por el método de Euler para h=0.001.]]
  
¿Cómo hay que elegir el paso de discretización temporal para que el método sea estable? 0.03 es estable
+
La estabilidad del método de Euler dependerá del comportamiento de la solución numérica del mismo cuando se perturba el valor de la condición inicial. Para comprobar cuándo el método es estable comparamos distintos resultados obtenidos al variar el valor de discretización (h) y observamos las diferencias. Se comprueba por tanto que el método es estable para h=0.05.
 +
[[Archivo:H05_leyenda_estable.png|600px|miniatura|centro|Solución para método estable para h=0.05.]]
  
 +
[[Archivo:H08_leyenda_inestable.png|600px|miniatura|centro|Solución método inestable para h=0.08.]]
 +
 +
=== Método del Trapecio ===
 +
 +
[[Archivo:Trap_estable_08_tn=3.png|450px|miniatura|derecha|Solución método del trapecio para tN=3.]]
 +
[[Archivo:Trap_estable_08_detalle_tN_0.2.png|450px|miniatura|derecha|Solución detalle método del trapecio con tN=0.2.]]
 +
 +
El método o regla del trapecio es un procedimiento basado en la integración numérica que permite calcular aproximadamente el valor de una integral definida. La regla se basa en aproximar el valor de la integral de \(f(x)\)  por el de la función lineal situado entre los puntos de inicio y final del intervalo de la función, aproximándolo al área de un trapecio. Utilizamos el siguiente código en Matlab para aplicar el método del trapecio (más exacto que Euler):
 +
 +
{{matlab|codigo=
 +
%Datos del problema:
 +
L=0.2;%Autoinductancia de la bobina
 +
R=5;  %Resistencia
 +
tau =L/R; %Definición
 +
E=20; %Valor de la fem
 +
t0=0; %Tiempo inicial
 +
tN=5*tau; %Tiempo de carga
 +
i0=0; %En t=0 está descargado
 +
funcion='100-i/0.04'; %i'=E/L-i/tau
 +
f=inline(funcion,'t','i');
 +
 +
%Discretización:
 +
h=0.08;
 +
N=round((tN-t0)/h);
 +
%Crear el vector t
 +
t=linspace(t0,tN,N+1);
 +
i=zeros(1,N+1);
 +
i(1)=i0;
 +
%Trapecio
 +
for k=1:N
 +
    i(k+1)=(1+h/(2*tau))\(i(k)+h/2*f(t(k),i(k))+h/2*E/L);
 +
end
 +
ir=E/R*(1-exp(-t./tau)); %solución analítica
 +
clf
 +
 +
%Representación:
 +
hold on
 +
plot(t,i,'b')
 +
plot(t,ir,'r')
 +
xlabel('tiempo')
 +
ylabel('intensidad')
 +
legend('Trapecio','Exacta','Location','Best');
 +
hold off
 +
}}
 +
 +
Al ser el método del Trapecio más efectivo que el de Euler, podemos utilizar un tamaño de paso de 0.08. Sin embargo, si utilizamos este tamaño para el de Euler el método es muy inestable, como hemos podido comprobar anteriormente.
 +
 +
Se ha elegido '''tN = 5*tau = 0.2s''' porque es el tiempo que tarda en cargarse la bobina, como puede apreciarse en la imagen siguiente:
 +
[[Archivo:Ppaa.jpg|300px|miniatura|centro|Solución por el método de trapecio.]]
 +
 +
Si aumentamos o disminuimos la constante del inductor:
 +
[[Archivo:Trap_distinto_L.png|800px|miniatura|centro|Intensidad para distintos valores de L.]]
 +
 +
== Interpretación de la ecuaciones en términos de las leyes de Kirchhoff ==
 +
 +
[[Image:Aaargggg.png|right|thumb|300px|Circuito RL]]
 +
Éste es el sistema de ecuaciones diferenciales corresponde al circuito de la figura de la derecha:
 +
<math>\left\{\begin{matrix}E(t)=R_1i_1(t)+L_2\frac{d}{dt}i_2(t)+R_2i_2(t)  [1]\\E(t)=R_1i_1(t)+L_1\frac{d}{dt}i_3(t)+R_3i_3(t)  [2]\\i_1(t)=i_2(t)+i_3(t)  [3]\end{matrix}\right.</math>
 +
 +
 +
La primera ley constituye la malla [2] a la cual consideramos la malla total. En ella se observa que la suma de potenciales (E(t)) se corresponde a la suma de intensidades por resistencias. Teniendo como resistencias R1 y R2 y L2 como inductancia. Y como intensidades i1, i2 e i3. La i2 corresponde a la intensidad que circula por la rama interior (1). La i3 circula, en su caso, por la rama restante al haber una diversificación de caminos, las intensidades son distintas.
 +
Estas intensidades están relacionadas con la primera ley de Kirchhoff, en la cual, la intensidad total o i1 corresponde a la suma de i2 e i3, según la ecuación [3] planteada en el nudo A.
 +
La ecuación [2] considera la malla 1. En ella la suma de potenciales es igual a E(t), la cual es equivalente a la suma de intensidades por resistencias. En este caso tenemos como resistencias R1 y R3 y como inductancia L3.
 +
 +
 +
<math>\left\{\begin{matrix}E(t)=R_1(i_2(t)+i_3(t))+L_2\frac{d}{dt}i_2(t)+R_2i_2(t)  [4]\\E(t)=R_1(i_2(t)+i_3(t))+L_1\frac{d}{dt}i_3(t)+R_3i_3(t)  [5]\end{matrix}\right.</math>
 +
 +
Agrupando y ordenando:
 +
 +
<math>\left\{\begin{matrix}E(t)=i_2(t)(R_1+R_2)+L_2\frac{d}{dt}i_2(t)+R_1i_3(t)  [6]\\E(t)=i_2(t)R_1+L_1\frac{d}{dt}i_3(t)+i_3(t)(R_1+R_3)  [7]\\i_2(0)=i_3(0)=0 [8]\end{matrix}\right.</math>
 +
 +
Las condiciones iniciales nos indican que el conector está abierto en el instante inicial. Al estar abierto no existe voltaje y, por lo tanto, tampoco corriente (circuito abierto).
 +
 +
== Resolución del sistema con datos ==
 +
 +
Resolvemos el sistema anterior mediante los siguientes valores: '''R1 = R2 = 6 Ω , R3 = 3 Ω , L1 = 0.3 H , L2 = 0.11 H , E(t) = 20 v constante.''' Al pasarlo a Matlab nos queda el programa siguiente utilizando los métodos de Euler explícito y trapecio implícito:
 +
 +
{{matlab|codigo=
 +
%Datos iniciales:
 +
to=0;
 +
tf=0.4;
 +
yo=[0,0]';
 +
N=1000;
 +
h=(tf-to)/N;
 +
t=to:h:tf;
 +
%f=@(y,t) [(y(1)*12+y(2)*6-20)/(-0.11);(y(1)*6+y(2)*9-20)/(-0.3)];
 +
f='[(y(1)*12+y(2)*6-20)/(-0.11);(y(1)*6+y(2)*9-20)/(-0.3)]';
 +
f=inline(f,'t','y');
 +
 +
%Euler:
 +
y= zeros(2,length(t));
 +
y(:,1)= yo;
 +
ytr=zeros(2,N+1);%Trapecio
 +
ytr(:,1)=yo;
 +
I=eye(2);%matriz identidad 2*2
 +
A=[-12/0.11, -6/0.11; -6/0.3, -9/0.3];%matriz de coeficientes
 +
T=[20/0.11;20/0.3];
 +
for i=1:N
 +
    y(:,i+1)=y(:,i)+h.*f(t(i),y(:,i));
 +
    ytr(:,1+i)=(I-h/2*A)\(ytr(:,i)+h/2*(f(t(i),y(:,i))+T));
 +
end
 +
i1=y(1,:)+y(2,:); %Hallamos i1
 +
i2=y(1,:); %extraemos i2 de la función y
 +
i3=y(2,:); %extraemos i3 de la función y
 +
y=[i1;y]; %Añadimos i1 a la función y
 +
 +
%Repetimos el proceso en el trapecio:
 +
i1tr=ytr(1,:)+ytr(2,:);
 +
i2tr=ytr(1,:); %extraemos i2 de la función y
 +
i3tr=ytr(2,:); %extraemos i3 de la función y
 +
ytr=[i1tr;ytr];  %Añadimos i1 a la función y trapecio
 +
hold on
 +
    title('Método de Euler');
 +
    plot(t,i1,'b')
 +
    plot(t,i2,'r')
 +
    plot(t,i3,'k')
 +
    legend('i1','i2','i3');
 +
 +
%Representamos:
 +
hold off
 +
figure
 +
hold on
 +
    title('Método del trapecio');
 +
    plot(t,i1tr,'.b')
 +
    plot(t,i2tr,'.r')
 +
    plot(t,i3tr,'.k')
 +
    legend('i1','i2','i3');
 +
hold off
 +
tt=0.4;
 +
p=round((tt-to)/h+1);
 +
[i1(p),i2(p)]
 +
[i1tr(p),i2tr(p)]
 +
}}
 +
 +
Comparamos los resultados de i1 e i2 para un diferencial de tiempo entre 0 y 0.4:
 +
[[Archivo:Intensidades.jpg|800px|miniatura|centro|Gráficos de intensidades con métodos de Euler y Trapecio.]]
 +
 +
i1 al ser la suma de i1 más i3 es mucho mayor que ambas llegando a valores de 2.5 A. i2 crece al principio más que i3 llegando a 1.25 A aproximadamente y luego desciende manteniéndose estable a menos de 1 A, mientras que i3 crece más lento pero nunca decrece y supera los 1.5 A.
 +
 +
Los valores en Euler y Trapecio no presentan grandes variaciones como podemos ver en las siguientes tablas. Los valores son cada vez más cercanos entre sí al aumentar el tiempo, llegando a diferenciarse en T=0.4 en 0.0001.
 +
[[Archivo:Tiemposint.png|800px|miniatura|centro|Gráficos comparativos de intensidades y tiempos.]]
 +
 +
== Modificación del Circuito RL ==
 +
 +
Puesto que el programa es idéntico al expuesto anteriormente no lo incluimos para evitar redundancias. La gráfica que obtenemos, en este caso al cambiar R3 de 3Ω a 9Ω, es la siguiente:
 +
[[Archivo:Intensidades_r3_9_apar_6.jpg|800px|miniatura|centro|Gráficos de intensidades con métodos de Euler y Trapecio.]]
 +
 +
Cambiando R3 = 3Ω por R3 = 9Ω se puede observar claramente el descenso en la intensidad i3 y el aumento de i2. La forma (pendiente) de la intensidad es la misma. La suma de i2 e i3, es decir, i1 disminuye de 2.5 A a poco mas de 2 A, el valor que alcanza cuando se vuelve constante. Esto es coherente con la ley de Ohm que dice que la resistencia R3 es inversamente proporcional a la corriente que circula sobre ella i3.
 +
 +
== Modificación de las intensidades ==
 +
Se modifican las condiciones iniciales de las intensidades en el sistema anterior de tal manera que i2 e i3 valen 1A en 0.02 s. Por lo tanto el valor inicial de las intensidades será:
 +
 +
{{matlab|codigo=
 +
to=0;
 +
tf=0.02;
 +
N=500;
 +
h=(tf-to)/N;
 +
yo=[1;1];
 +
%t=to:h:tf;
 +
t=linspace(to,tf,N+1);
 +
y=zeros(2,N+1);
 +
y(:,N+1)= yo;
 +
A=[-12/0.11, -6/0.11; -6/0.3, -9/0.3];
 +
T=[20/0.11;20/0.3];
 +
f=inline('[(y(1)*12+y(2)*6-20)/(-0.11);(y(1)*6+y(2)*15-20)/(-0.3)]','t','y');
 +
for i=N+1:-1:2
 +
y(:,i-1)=y(:,i)-h.*f(t(i),y(:,i));
 +
end
 +
y(:,1)
 +
plot(t,y,'.')
 +
legend('i2(t)','i3(t)')
 +
}}
  
 +
[[Archivo:Intensidades_en_0.jpg|600px|miniatura|centro|Intensidades a partir del cambio de condiciones iniciales.]]
  
 +
Los valores iniciales de i2 e i3 serán -0.2935 y 0.8893 respectivamente, para que se cumpla que i2 e i3 en 0.02 sea 1 A. Los signos negativos significan que la corriente ha sufrido un cambio de sentido.
  
  

Revisión actual del 14:10 27 abr 2016


1 Introducción

Un circuito RL es un circuito eléctrico que contiene una resistencia y una bobina en serie, además de una fuente de alimentación. Se dice que la bobina se opone transitoriamente al establecimiento de una corriente en el circuito. En una resistencia R, la Ley de Ohm establece que:

[math]i(t)={v(t)\over R}[/math] siendo [math]i(t)[/math] la intensidad de corriente en amperios ([math]A[/math]), [math]v(t)[/math] el voltaje en voltios ([math]V[/math]) y [math]R[/math] el coeficiente de resistencia en ohmios ([math]Ω[/math]).

En un inductor L, la Ley de Faraday impone:

[math]v(t)=L {d\over dt} i(t)[/math] donde [math]L[/math] es el coeficiente de autoinducción en henrios ([math]H[/math]).

También tenemos en cuenta las Leyes de Kirchhoff, que establecen el comportamiento de los circuitos:

  1. Ley de corrientes: En cada nodo, la suma de corrientes que entra es igual a la que sale.
  2. Ley de tensiones: En cada ciclo cerrado o malla, la suma de diferencias de potencial es nula.


Circuitos Eléctricos RL

2 Ley de Kirchhoff de Voltaje o Tensiones

La ley de voltaje de Kirchhoff indica que la suma de voltajes alrededor de una trayectoria o circuito cerrado debe ser cero. \[{d\over dt}i(t)+{R\over L}i(t)-{E(t)\over L}=0\]

centro

Por estar en serie y aplicando esta ley, podemos ver que la tensión total es la tensión en la resistencia (R) más la tensión en la bobina (L). Por lo que la f.e.m. es igual al voltaje de la bobina más la resistencia.

centro

Aplicando la ley de Ohm a la resistencia y la de Faraday a la bobina, se obtiene la siguiente ecuación:

centro

2.1 Cálculo analítico y representación

Gráfico representación analítica intensidad-tiempo.

Suponiendo que en el instante t=0 el circuito pasa de estar abierto a cerrado, obtenemos la intensidad para cada instante t>0 teniendo en cuenta que la fuente de alimentación tiene un voltaje constante de E=20V, L=0.2H y R=5Ω. El cálculo analítico de la intensidad quedará resuelto de la siguiente manera, mediante un Problema de Valor Inicial (P.V.I.) o de Cauchy:


Cálculos analíticos.

2.2 Método de Euler

Aplicamos el método de Euler basado en la aproximación del valor de la función a la recta tangente en cada punto conocido, mediante un problema de valor inicial (P.V.I.) y lo comparamos con los resultados que hemos obtenido analíticamente.

Gráfico intensidad respecto del tiempo cuando el circuito pasa de estar abierto a cerrado aplicando el método numérico.
%Datos del problema: 
L=0.2;%Autoinductancia de la bobina 
R=5;%Resistencia
tau =L/R; %Definición
E=20; %Valor de la fem
t0=0;%Tiempo
tN=5*tau; %Tiempo de carga
i0=0; %En t=0 está descargado
funcion='100-i/0.04'; %i'=E/L-i/tau
f=inline(funcion,'t','i');

%Discretización:
h=0.001;
N=round((tN-t0)/h);
%Crear el vector t
t=linspace(t0,tN,N+1);
i=zeros(1,N+1);
i(1)=i0;
%Euler
for k=1:N
    i(k+1)=i(k)+h*f(t(k),i(k));
end
ir=E/R*(1-exp(-t./tau)); %solción analítica
clf

%Representación:
hold on
plot(t,i,'b')
plot(t,ir,'r')
xlabel('tiempo')
ylabel('intensidad')
legend('Euler','Exacta','Location','Best');
hold off


En la gráfica podemos ver que efectivamente el tiempo de carga es 5*tau.

Solución por el método de Euler para h=0.01.
Solución por el método de Euler para h=0.001.

La estabilidad del método de Euler dependerá del comportamiento de la solución numérica del mismo cuando se perturba el valor de la condición inicial. Para comprobar cuándo el método es estable comparamos distintos resultados obtenidos al variar el valor de discretización (h) y observamos las diferencias. Se comprueba por tanto que el método es estable para h=0.05.

Solución para método estable para h=0.05.
Solución método inestable para h=0.08.

2.3 Método del Trapecio

Solución método del trapecio para tN=3.
Solución detalle método del trapecio con tN=0.2.

El método o regla del trapecio es un procedimiento basado en la integración numérica que permite calcular aproximadamente el valor de una integral definida. La regla se basa en aproximar el valor de la integral de \(f(x)\) por el de la función lineal situado entre los puntos de inicio y final del intervalo de la función, aproximándolo al área de un trapecio. Utilizamos el siguiente código en Matlab para aplicar el método del trapecio (más exacto que Euler):

%Datos del problema:
L=0.2;%Autoinductancia de la bobina 
R=5;  %Resistencia
tau =L/R; %Definición
E=20; %Valor de la fem
t0=0; %Tiempo inicial
tN=5*tau; %Tiempo de carga
i0=0; %En t=0 está descargado
funcion='100-i/0.04'; %i'=E/L-i/tau
f=inline(funcion,'t','i');

%Discretización:
h=0.08;
N=round((tN-t0)/h);
%Crear el vector t
t=linspace(t0,tN,N+1);
i=zeros(1,N+1);
i(1)=i0;
%Trapecio
for k=1:N
    i(k+1)=(1+h/(2*tau))\(i(k)+h/2*f(t(k),i(k))+h/2*E/L);
end
ir=E/R*(1-exp(-t./tau)); %solución analítica
clf

%Representación:
hold on
plot(t,i,'b')
plot(t,ir,'r')
xlabel('tiempo')
ylabel('intensidad')
legend('Trapecio','Exacta','Location','Best');
hold off


Al ser el método del Trapecio más efectivo que el de Euler, podemos utilizar un tamaño de paso de 0.08. Sin embargo, si utilizamos este tamaño para el de Euler el método es muy inestable, como hemos podido comprobar anteriormente.

Se ha elegido tN = 5*tau = 0.2s porque es el tiempo que tarda en cargarse la bobina, como puede apreciarse en la imagen siguiente:

Solución por el método de trapecio.

Si aumentamos o disminuimos la constante del inductor:

Intensidad para distintos valores de L.

3 Interpretación de la ecuaciones en términos de las leyes de Kirchhoff

Circuito RL

Éste es el sistema de ecuaciones diferenciales corresponde al circuito de la figura de la derecha: [math]\left\{\begin{matrix}E(t)=R_1i_1(t)+L_2\frac{d}{dt}i_2(t)+R_2i_2(t) [1]\\E(t)=R_1i_1(t)+L_1\frac{d}{dt}i_3(t)+R_3i_3(t) [2]\\i_1(t)=i_2(t)+i_3(t) [3]\end{matrix}\right.[/math]


La primera ley constituye la malla [2] a la cual consideramos la malla total. En ella se observa que la suma de potenciales (E(t)) se corresponde a la suma de intensidades por resistencias. Teniendo como resistencias R1 y R2 y L2 como inductancia. Y como intensidades i1, i2 e i3. La i2 corresponde a la intensidad que circula por la rama interior (1). La i3 circula, en su caso, por la rama restante al haber una diversificación de caminos, las intensidades son distintas. Estas intensidades están relacionadas con la primera ley de Kirchhoff, en la cual, la intensidad total o i1 corresponde a la suma de i2 e i3, según la ecuación [3] planteada en el nudo A. La ecuación [2] considera la malla 1. En ella la suma de potenciales es igual a E(t), la cual es equivalente a la suma de intensidades por resistencias. En este caso tenemos como resistencias R1 y R3 y como inductancia L3.


[math]\left\{\begin{matrix}E(t)=R_1(i_2(t)+i_3(t))+L_2\frac{d}{dt}i_2(t)+R_2i_2(t) [4]\\E(t)=R_1(i_2(t)+i_3(t))+L_1\frac{d}{dt}i_3(t)+R_3i_3(t) [5]\end{matrix}\right.[/math]

Agrupando y ordenando:

[math]\left\{\begin{matrix}E(t)=i_2(t)(R_1+R_2)+L_2\frac{d}{dt}i_2(t)+R_1i_3(t) [6]\\E(t)=i_2(t)R_1+L_1\frac{d}{dt}i_3(t)+i_3(t)(R_1+R_3) [7]\\i_2(0)=i_3(0)=0 [8]\end{matrix}\right.[/math]

Las condiciones iniciales nos indican que el conector está abierto en el instante inicial. Al estar abierto no existe voltaje y, por lo tanto, tampoco corriente (circuito abierto).

4 Resolución del sistema con datos

Resolvemos el sistema anterior mediante los siguientes valores: R1 = R2 = 6 Ω , R3 = 3 Ω , L1 = 0.3 H , L2 = 0.11 H , E(t) = 20 v constante. Al pasarlo a Matlab nos queda el programa siguiente utilizando los métodos de Euler explícito y trapecio implícito:

%Datos iniciales:
to=0;
tf=0.4;
yo=[0,0]';
N=1000;
h=(tf-to)/N;
t=to:h:tf;
%f=@(y,t) [(y(1)*12+y(2)*6-20)/(-0.11);(y(1)*6+y(2)*9-20)/(-0.3)];
f='[(y(1)*12+y(2)*6-20)/(-0.11);(y(1)*6+y(2)*9-20)/(-0.3)]';
f=inline(f,'t','y');

%Euler:
y= zeros(2,length(t));
y(:,1)= yo;
ytr=zeros(2,N+1);%Trapecio
ytr(:,1)=yo;
I=eye(2);%matriz identidad 2*2
A=[-12/0.11, -6/0.11; -6/0.3, -9/0.3];%matriz de coeficientes
T=[20/0.11;20/0.3];
for i=1:N
    y(:,i+1)=y(:,i)+h.*f(t(i),y(:,i));
    ytr(:,1+i)=(I-h/2*A)\(ytr(:,i)+h/2*(f(t(i),y(:,i))+T));
end
i1=y(1,:)+y(2,:); %Hallamos i1
i2=y(1,:); %extraemos i2 de la función y
i3=y(2,:); %extraemos i3 de la función y
y=[i1;y]; %Añadimos i1 a la función y

%Repetimos el proceso en el trapecio:
i1tr=ytr(1,:)+ytr(2,:);
i2tr=ytr(1,:); %extraemos i2 de la función y
i3tr=ytr(2,:); %extraemos i3 de la función y
ytr=[i1tr;ytr];   %Añadimos i1 a la función y trapecio
hold on
    title('Método de Euler');
    plot(t,i1,'b')
    plot(t,i2,'r')
    plot(t,i3,'k')
    legend('i1','i2','i3');

%Representamos:
hold off
figure
hold on
    title('Método del trapecio');
    plot(t,i1tr,'.b')
    plot(t,i2tr,'.r')
    plot(t,i3tr,'.k')
    legend('i1','i2','i3');
hold off
tt=0.4;
p=round((tt-to)/h+1);
[i1(p),i2(p)]
[i1tr(p),i2tr(p)]


Comparamos los resultados de i1 e i2 para un diferencial de tiempo entre 0 y 0.4:

Gráficos de intensidades con métodos de Euler y Trapecio.

i1 al ser la suma de i1 más i3 es mucho mayor que ambas llegando a valores de 2.5 A. i2 crece al principio más que i3 llegando a 1.25 A aproximadamente y luego desciende manteniéndose estable a menos de 1 A, mientras que i3 crece más lento pero nunca decrece y supera los 1.5 A.

Los valores en Euler y Trapecio no presentan grandes variaciones como podemos ver en las siguientes tablas. Los valores son cada vez más cercanos entre sí al aumentar el tiempo, llegando a diferenciarse en T=0.4 en 0.0001.

Gráficos comparativos de intensidades y tiempos.

5 Modificación del Circuito RL

Puesto que el programa es idéntico al expuesto anteriormente no lo incluimos para evitar redundancias. La gráfica que obtenemos, en este caso al cambiar R3 de 3Ω a 9Ω, es la siguiente:

Gráficos de intensidades con métodos de Euler y Trapecio.

Cambiando R3 = 3Ω por R3 = 9Ω se puede observar claramente el descenso en la intensidad i3 y el aumento de i2. La forma (pendiente) de la intensidad es la misma. La suma de i2 e i3, es decir, i1 disminuye de 2.5 A a poco mas de 2 A, el valor que alcanza cuando se vuelve constante. Esto es coherente con la ley de Ohm que dice que la resistencia R3 es inversamente proporcional a la corriente que circula sobre ella i3.

6 Modificación de las intensidades

Se modifican las condiciones iniciales de las intensidades en el sistema anterior de tal manera que i2 e i3 valen 1A en 0.02 s. Por lo tanto el valor inicial de las intensidades será:

to=0; 
tf=0.02;
N=500;
h=(tf-to)/N;
yo=[1;1];
%t=to:h:tf;
t=linspace(to,tf,N+1);
y=zeros(2,N+1);
y(:,N+1)= yo;
A=[-12/0.11, -6/0.11; -6/0.3, -9/0.3];
T=[20/0.11;20/0.3];
f=inline('[(y(1)*12+y(2)*6-20)/(-0.11);(y(1)*6+y(2)*15-20)/(-0.3)]','t','y');
for i=N+1:-1:2
 y(:,i-1)=y(:,i)-h.*f(t(i),y(:,i));
end
y(:,1)
plot(t,y,'.')
legend('i2(t)','i3(t)')


Intensidades a partir del cambio de condiciones iniciales.

Los valores iniciales de i2 e i3 serán -0.2935 y 0.8893 respectivamente, para que se cumpla que i2 e i3 en 0.02 sea 1 A. Los signos negativos significan que la corriente ha sufrido un cambio de sentido.