Diferencia entre revisiones de «Circuitos RL (grupo B)»

Revisión actual del 00:37 7 mar 2014

El circuito eléctrico mas simple esta compuesto de una resistencia, un inductor o bobina y una fuente de alimentación.

• En una resistencia R, la Ley de Ohm establece:

 [math]i(t)={V(t)\over R}[/math]

• En un inductor L la Ley de Faraday dice:

[math] V(t)=L\cdot i'(t)[/math]

Las leyes de Kirchhoff establecen el comportamiento de los circuitos:

1. Ley de corriente: en cada nudo la suma de corrientes que entra coincide con la suma de corrientes que sale.
2. Ley de tensiones: en cada ciclo cerrado, la suma de las diferencias de potencial eléctrico es nula.

1 Segunda ley de Kirchhoff para un circuito RL simple.

1.1 Ecuación diferencial-Resolución y representación analítica.

Basándonos en la segunda ley de Kirchhoff sobre el circuito que observamos en la imagen, podemos trabajar con la siguiente ecuación diferencial.

Esta vez vamos a obtener la solución del valor concreto de un punto de una función hallada a partir de una ecuación diferencial. Será así de nuevo el caso de una ecuación diferencial que define las relaciones entre las corrientes y sus derivadas en un circuito RL de varias resistencias y bobinas. Volveremos a utilizar el segundo circuito mostrado en este artículo:

Circuito RL simple

Al suponer que en el instante t=0 el circuito está abierto, podemos establecer para el problema de Cauchy el valor inicial [math] i_0(t)=0 [/math]. Considerando que la alimentación posee un voltaje constante V(t)=20V, L=0.2 y R=5Ω, llegamos al problema:

                             [math]\left\{\begin{matrix}\ i'(t)+{25}i(t)-{100}=0 \\i(0)=0\end{matrix}\right.[/math]

La resolución de éste mediante [math]i \cdot e^{\int{\frac{R}{L}dt}} = \int{e^{\frac{R \cdot t}{L}}dt} \hspace{0.5cm} \rightarrow \hspace{0.5cm} [/math] [math]i \cdot e^{\frac{R \cdot t}{L}} = \frac{V}{R}e^{\frac{R \cdot t}{L}}+C[/math]
nos da:

[math] i(t)=4-4e^{-25t} [/math]

La gráfica de la función nos muestra la evolución de la intensidad. En ella se puede apreciar una rápida subida en los primeros instantes que se ralentiza hasta llegar a ser prácticamente constante con un valor de 4. Por lo que podemos afirmar que la intensidad se estabiliza en un corto periodo de tiempo.

Gráfica de la intensidad en un circuito RL simple

fplot('4-4*exp(-25*t)',[0,0.5,0,5]);
xlabel('Tiempo')
ylabel('Intensidad(t)')

1.2 Método de Euler.

Este método se basa en la posibilidad de aproximar una función original mediante rectas tangentes partiendo de un punto dado. Sabiendo que la recta tangente de una función en un punto es la derivada en dicho punto, este método consiste en calcular derivadas de la función en diferentes puntos pertenecientes a ella. Estos puntos serán tomados con la misma distancia a lo largo de la recta, a esta distancia la denominamos paso de discretización.

Aproximación de una función por el método de Euler

A continuación, vamos a calcular la función por el método de Euler, resolviendo la ecuación con diferentes pasos. De esta forma, podremos observar y analizar los distintos resultados obtenidos al modificarlos, hasta llegar a una conclusión evidente. Además calcularemos el error cometido por la aproximación del método para ser más precisos en nuestras observaciones.

Comenzamos a analizar este método elaborando un programa con un paso de 0'01, que podemos visualizar a continuación.

%Definimos las condiciones iniciales.
clear all
t0=0;
tN=1;

%Establecemos un paso h=0'01.
h=0.01;
t=t0:h:tN;
N=(tN-t0)/h;

%Resolvemos la ecuación.
ii=0;
i(1)=ii;
for n=1:N;
  ii=ii+h*(100-25*ii);
  i(n+1)=ii;
end

%Dibujamos la gráfica.
figure(1)
plot(t,i,'-g','linewidth',2);
xlabel('Tiempo');
ylabel('Intensidad(t)');
axis([-0.1,1,0,5])

%Calculo error y lo dibujamos.
r=4-4*exp(-25*t);
error=abs(r-i);
figure(2)
plot(t,error,'-r','linewidth',2)
xlabel('Tiempo');
ylabel('Error');
axis([-0.1,1,-0.05,0.25])

Un paso de discretización de 0'01 proporciona un error muy pequeño, lo que aporta una aproximación de la función muy similar a la real, como podemos observar en la primera figura. Analizando el error que nos ofrece dicha aproximación (imagen de la derecha),su valor no supera el 0'2, por tanto, podemos afirmar analíticamente la exactitud de la aproximación. Además, podemos observar que el error se anula al llegar a t=0'25, momento en el que coincide la función aproximada por Euler con la función real.

Resultado de la intensidad del circuito RL simple por el método de Euler para h=0'01

Gráfica del error cometido en el cálculo de la intensidad por el método de Euler para h=0'01

Una vez analizado el resultado anterior, resolvemos el problema con un paso mayor.

%Resolucion por Euler
%Definimos las condiciones iniciales.
clear all
t0=0;
tN=1;

%Probamos con un paso de discretización mayor, h=0'05.
h=0.05;
t=t0:h:tN;
N=(tN-t0)/h;

%Resolvemos la ecuación.
ii=0;
i(1)=ii;
for n=1:N;
  ii=ii+h*(100-25*ii);
  i(n+1)=ii;
end

%Creamos el dibujo.
figure(1)
plot(t,i,'-g','linewidth',2);
xlabel('Tiempo');
ylabel('Intensidad(t)');
axis([-0.1,1,0,5.1])

%Calculo error y dibujamos.
r=4-4*exp(-25*t);
error=abs(r-i);
figure(2)
plot(t,error,'-r','linewidth',2)
axis([-0.1,1,-0.05,2.25])
xlabel('Tiempo');
ylabel('Error');

Podemos observar que al aumentar el paso de discretización a 0'05 la aproximación de la curva se ha convertido en una combinación de rectas que asemejan la evolución de la intensidad, ya que sigue manteniendo una rápida subida inicial. A pesar de ello, el error ha aumentado su valor en un 1000% en el periodo transitorio, pero vuelve a anularse en el permanente ya que a partir del punto de estabilización los valores aproximados coincide con los reales.El cambio producido en el incremento inicial de la intensidad se debe a la brevedad de éste, por lo que resulta poco preciso trabajar con pasos mayores.

Resultado de la intensidad del circuito RL simple por el método de Euler para h=0'05

Gráfica del error cometido en el cálculo de la intensidad por el método de Euler para h=0'05

Para estar más seguro de la afirmación hecha hace unos instantes, en la que hemos asegurado que es poco preciso trabajar con pasos mayores, vamos a resolver y comparar en la misma gráfica el problema resuelto con distintos valores para la h, el paso.

%Comparamos las gráficas con diferentes pasos y la función real.
clear all

%Guardamos el valor para cada paso.
t0=0;
tN=1;
h1=0.05;
t1=t0:h1:tN;
N1=(tN-t0)/h1;
h2=0.01;
t2=t0:h2:tN;
N2=(tN-t0)/h2;
h3=0.075;
t3=t0:h3:tN;
N3=(tN-t0)/h3;

%Creamos un bucle para cada h y resolver las cuatro ecuaciones.
ii=0;
i(1)=ii;
for n=1:N1;
  ii=ii+h1*(100-25*ii);
  i(n+1)=ii;
end
yy=0;
y(1)=yy;
for m=1:N2;
  yy=yy+h2*(100-25*yy);
  y(m+1)=yy;
end
jj=0;
j(1)=jj;
for m=1:N3;
  jj=jj+h3*(100-25*jj);
  j(m+1)=jj;
end

%Hacemos el dibujo con todas las gráficas a la vez para hacer una comparación global de estas y la función real.
hold on 
fplot('4-4*exp(-25*t2)',[0,1,0,5],'r');
plot(t1,i,'-b','linewidth',1);
plot(t2,y,'-y','linewidth',1);
plot(t3,j,'-g','linewidth',1);
legend('Función real','h=0.05', 'h=0.01', 'h=0.075');
title('Comparativa');
xlabel('Tiempo');
ylabel('Intensidad(t)');
axis([-0.1,1,-0.05,8])
hold off

%Quitamos el paso mayor para ver con mayor exactitud como se difieren las aproximaciones más correctas.
figure (2)
hold on 
fplot('4-4*exp(-25*t2)',[0,1,0,5],'r');
plot(t1,i,'-b','linewidth',1);
plot(t2,y,'-g','linewidth',1);
legend('Función real','h=0.05', 'h=0.01');
title('Comparativa');
xlabel('Tiempo');
ylabel('Intensidad(t)');
axis([0,0.5,-0.05,5])
hold off

Al ensayar con diferentes pasos de discretización vemos como la aproximación de la gráfica se va degenerando, hasta el punto en el que si se sobrepasa 0'05, la volatilidad de las oscilaciones llega a ser tal que la intensidad tarda un periodo sensiblemente mayor en estabilizarse y el error no se anula, ya que los valores son muy lejanos a los reales. Como la subida de intensidad se produce muy rápidamente, una h≥0'05 es, como podemos ver la imagen de la izquierda, muy poco precisa. Si se mantiene en el intervalo [0'02,0'05] se elimina la forma exponencial de la función (aparecen en su lugar fluctuaciones rectilíneas) pero se sigue reflejando la evolución de la intensidad hasta que se vuelve constante. Con un paso de 0'01, se consigue una gráfica exponencial con valores ligeramente superiores a los reales, por tanto podemos afirmar cuanto más pequeño sea el paso, más exacta y precisa será la aproximación ya que se trata de un aumento prácticamente instantáneo.

Comparativa de la función real con aproximaciones con diferentes pasos de discretización.

Comparativa de la función real con aproximaciones con pasos de discretización más precisos.

1.3 Método del Trapecio.

El método del trapecio se basa en aproximar una función integrando. Para ello utilizamos un número N de trapecios (en la imagen inferior hemos aproximado mediante 4 trapecios). Las integrales a resolver para poder poder aproximar una función por este método tendrán la distancia de sus límites de integración constantes para todas las integrales, a esta distancia entre los límites de integración lo llamamos paso. En el desarrollo de este método se supone que la función es continua en el intervalo en que se quiere aproximar.

Función Y(x), el área bajo la curva, se puede aproximar mediante n trapecios

A continuación, vamos a resolver el problema mediante el método del trapecio con un paso de 0'05.

%Resolucion por Trapecio.
%Definimos las condiciones iniciales.
t0=0;
tN=1;

%Al tratarse del método de trapecio podemos usar un paso más grande.
h=0.05;
t=t0:h:tN;
N=(tN-t0)/h;
i=zeros(1,N+1);
i(1)=0;

%Resolvemos la ecuación.
for n=1:N;
    i(n+1)=((200*h)/(2+25*h))+(i(n)*((2-25*h)/(2+25*h)));
end

%Dibujamos.
figure (1)
plot(t,i,'-r','linewidth',2)
xlabel('Tiempo')
ylabel('Intensidad(t)')
axis([0,1,-0.05,4.25])

%Cambiamos los ejes del dibujo y añadimos líneas auxiliares para poder ver con mayor 
%precisión el punto en el que la intensidad se hace constante.
figure (2)
plot(t,i,'-b','linewidth',2)
xlabel('Tiempo')
ylabel('Intensidad(t)')
axis([0,0.4,2,4.1])
grid on

%Dibujamos el error
r=4-4*exp(-25*t);
error=abs(i-r);
figure(3)
plot(t,error,'-g','linewidth',2)
xlabel('Tiempo')
ylabel('Intensidad(t)')

Resultado de la intensidad del circuito RL simple por el método del trapecio.

Visualización del punto en que se hace constante la intensidad.

Vemos en la imagen inferior que el error cometido en el método del Trapecio para h=0'05 es similar al cometido con el método de Euler y un paso menor. Por tanto podemos afirmar que éste es más preciso para pasos mayores.

Visualización del punto en que se hace constante la intensidad.

La representación de la función por el método del Trapecio nos permite utilizar un paso de discretización de 0'5 quedando como resultado una aproximación muy similar a la original (podemos observar ésto en la imagen de la izquierda). Una ampliación de dicha aproximación nos permite ver con exactitud el punto en el cual la intensidad se hace constante,(0'25,4), es decir t=0'25 e i(t)=4.

Además, si comparamos las gráficas obtenidas al aproximar la función por el método de Euler y del Trapecio con el mismo paso podemos afirmar que el método del Trapecio hace una aproximación más precisa que el método de Euler.

Método del Trapecio (h=0.05).

Método de Euler (h=0'05)

1.3.1 Efectos en la función al variar la constante del inductor.

La bobina es un elemento pasivo del circuito capaz de almacenar energía magnética.

Una vez explicado la función del inductor (bobina) en el circuito, nos disponemos a estudiar las alteraciones producidas en la función cuando se disminuye o aumenta la constante del inductor, L.

A continuación se muestran las gráficas resultantes al modificar dicho valor.

%Resolucion por Trapecio.
%Definimos las condiciones iniciales.
t0=0;
tN=1;
h=0.05;
t=t0:h:tN;
N=(tN-t0)/h;
i=zeros(N+1);
i(1)=0;

%Cambiamos la L y resolvemos su ecuación correspondiente.
L=3.2;
for n=1:N;
    i(n+1)=((40*h)/(2*L+5*h))+i(n)*((2*L-5*h)/(2*L+5*h));
end
figure (1)
plot(t,i,'-b','linewidth',2)
xlabel('Tiempo')
ylabel('Intensidad(t)')
title('L=3.2')
grid on

%Cambiamos la L y resolvemos su ecuación correspondiente.
y=zeros(N+1);
y(1)=0;
L=0.8;
for m=1:N;
    y(m+1)=((40*h)/(2*L+5*h))+y(m)*((2*L-5*h)/(2*L+5*h));
end
figure (2)
plot(t,y,'-r','linewidth',2)
xlabel('Tiempo')
ylabel('Intensidad(t)')
title('L=0.8')
grid on

%Cambiamos la L y resolvemos su ecuación correspondiente.
j=zeros(N+1);
j(1)=0;
L=0.05;
for o=1:N;
    j(o+1)=((40*h)/(2*L+5*h))+j(o)*((2*L-5*h)/(2*L+5*h));
end
figure (3)
plot(t,j,'-g','linewidth',2)
xlabel('Tiempo')
ylabel('Intensidad(t)')
title('L=0.05')
grid on

%Cambiamos la L y resolvemos su ecuación correspondiente.
k=zeros(N+1);
k(1)=0;
L=0.0125;
for l=1:N;
    k(l+1)=((40*h)/(2*L+5*h))+k(l)*((2*L-5*h)/(2*L+5*h));
end
figure (4)
plot(t,k,'-c','linewidth',2)
xlabel('Tiempo')
ylabel('Intensidad(t)')
title('L=0.0125')
grid on

• Como podemos observar a medida que incrementamos el valor de L, se acentúa la 'exponencialidad' de la función hasta llegar a aproximarse a una recta. A medida que aumenta el parecido a una recta, el tiempo de estabilización se extiende por ser la inductancia mayor, aunque el valor de la intensidad cuando es constante se mantiene.
• De igual forma al reducirlo, la subida de la intensidad hasta llegar a una cifra constante es mucho mayor, aumenta el periodo transitorio, pero de forma más moderada que en los casos anteriores y manteniendo la constancia en i(t)=4; sin embargo, la función deja de ser exponencial convirtiéndose en fluctuaciones de rectas que van reduciendo su volatilidad con el paso del tiempo.

L=3'2.

L=0'8.

L=0'0125.

L=0'5.

2 Leyes de Kirchhoff aplicada a un circuito RL de dos mallas.

Tras estudiar el circuito RL más simple (de una sola malla), ahora estudiaremos un circuito de dos mallas, por lo que como antes nos aparecía una sola ecuación diferencial, ahora nos aparecerá un sistema de dos ecuaciones diferenciales (una por cada malla).

Circuito RL de varias resistencias y bobinas

Según las leyes de Kirchhoff, el sistema de ecuaciones diferenciales correspondiente al circuito 2 es el siguiente:

[math]\left\{\begin{matrix}\ [1] E(t) = R_1\cdot i_1(t) + L_2\cdot \frac{d}{dt}\cdot i_2(t) + R_2\cdot i_2(t)\\ [2] E(t) = R_1\cdot i_1(t) + L_1\cdot \frac{d}{dt}\cdot i_3(t) + R_3\cdot i_3(t)\\ [3] i_1(t) = i_2(t) + i_3(t)\end{matrix}\right.[/math]

[1]Correspondiente al recorrido exterior del circuito.
[2]Correspondiente a la malla izquierda, (malla 1).
[3]Correspondiente al nudo de la parte superior del circuito(Ley de corrientes de Kirchhoff).

Sustituyendo la tercera ecuación en las otras dos, se obtiene el sistema en términos de i2(t) e i3(t) matricialmente::[math]\frac{d}{dt}\cdot \begin{pmatrix} i_2\\i_3 \end{pmatrix} = - \begin{pmatrix} \frac{R_1+R_2}{L_2}&\frac{R_1}{L_2} \\ \frac{R_1}{L_1}&\frac{R_1+R_3}{L_1}\end{pmatrix} \cdot \begin{pmatrix} i_2\\i_3 \end{pmatrix} + \begin{pmatrix} \frac{E}{L_2} \\ \frac{E}{L_1} \end{pmatrix}[/math]
A partir de las condiciones iniciales i2(0)=i3(0)=0 se puede interpretar que la corriente en el circuito es nula en el momento en el cual se conecta el generador.

2.1 Cálculo de las intensidades de un circuito con varias resistencias y bobinas a partir del Método de Euler y del Trapecio.

Ya se ha visto que la expresión de las intensidades i2 e i3 de este circuito RL vienen dadas por dos ecuaciones diferenciales que conforman un sistema. Este sistema de ecuaciones diferenciales se puede resolver numéricamente a partir del método de Euler y del trapecio, los cuales se basan en la aproximación de las áreas de esa ecuación diferencial al valor real de la función en cada momento.

Partiendo de unos datos iniciales iguales a:

[math] R_1=R_2=6Ω [/math] , [math] R_3=3Ω [/math] , [math] L_1=0,3Hz [/math] , [math] L_2=0,11Hz [/math] y [math] E(t)=20 V [/math]

:[math]\frac{d}{dt}\cdot \begin{pmatrix} i_2\\i_3 \end{pmatrix} = - \begin{pmatrix} \frac{6+6}{0,11}&\frac{6}{0,11} \\ \frac{6}{0,3}&\frac{6+3}{0,3}\end{pmatrix} \cdot \begin{pmatrix} i_2\\i_3 \end{pmatrix} + \begin{pmatrix} \frac{20}{0,11} \\ \frac{20}{0,3} \end{pmatrix}[/math]

Resolvemos las ecuaciones diferenciales:

• MÉTODO DE EULER:

OCTAVE

DATOS:

clc;clear all;
    %t0=tiempo inicial
t0=0;
    %tN=tiempo máximo a evaluar
tN=0.04;
    %soluciones en y0 para las variables i2 e i3
y0=0;
    %N=intervalos en los que se dividen
N=25;
    %h=salto de intervalo en intervalo
h=(tN-t0)/N;
    %i1,i2,i3=vectores solución desde t0 a tN de h en h
i1=t0:h:tN;

i2=t0:h:tN;

i3=t0:h:tN;

    %y=vector de ceros y que va a almacenar las soluciones.Tendrá N+1 términos porque el primero es el cero y el ultimo el N
    %la intensidad dos es y(:,1),es decir, la fila 1.
    %la intensidad tres es y(:,2),es decir, la fila 2.
    %la intensidad uno es y(:,3),es decir, la fila 3.
y=zeros(N+1,3);
y(1,1)=y0;
y(1,2)=y0;
y(1,3)=y0;
%AHORA SE EMPIEZA A RESOLVER POR EULER EL SISTEMA.
%y(n+1)=y(n)+h*f(t(n),y(n))
%f(t,i2)=((E/L2)-R1*(i2+i3)+R2*i2)/L2
%f(t,i3)=((E/L1)-R1*(i2+i3)+R3*i3)/L1
for n=1:N
    y(n+1,1)=y(n,1)+h*((20/0.11)-(6*(y(n,1)+y(n,2))+6*y(n,1))/0.11);
    y(n+1,2)=y(n,2)+h*((20/0.3)-(6*(y(n,1)+y(n,2))+3*y(n,2))/0.3);
    y(n+1,3)=y(n+1,1)+y(n+1,2);
end
hold on
plot(i1,y(:,3),'+')
plot(i2,y(:,1),'r+')
plot(i3,y(:,2),'g+')
xlabel('Tiempo');
ylabel('Intensidad(t)');
legend('i1', 'i2', 'i3');
grid on

• MÉTODO DEL TRAPECIO:

OCTAVE

DATOS:
clc;clear all;
    %tiempo inicial
t0=0;
    %tiempo máximo a evaluar
tN=0.04;
    %soluciones en y0 para las variables i2 e i3
y0=0;
    %intervalos en los que se dividen
N=25;
	%salto de intervalo en intervalo
h=(tN-t0)/N;
    %vectores solución desde t0 a tN de h en h
i1=t0:h:tN;
i2=t0:h:tN;
i3=t0:h:tN;
    %vector de ceros y que va a almacenar las soluciones
    %tendra N+1 términos porque el primero es el cero y el ultimo el N
    %la intensidad dos es y(:,1)
    %la intensidad tres es y(:,2)
    %la intensidad uno es y(:,3)
y=zeros(N+1,3);
y(1,1)=y0;
y(1,2)=y0;
y(1,3)=y0;
PARA RESOLVER POR EL MÉTODO DEL TRAPECIO:
%y(n+1)=y(n)+(h/2)*(f(t(n),y(n)+f(t(n+1),y(n+1))
%f(t,i2)=((E/L2)-R1*(i2+i3)+R2*i2)/L2
%f(t,i3)=((E/L1)-R1*(i2+i3)+R3*i3)/L1
for n=1:N
    y(n+1,1)=y(n,1)+(h/2)*(((20/0.11)-(6*(y(n,1)+y(n,2))+6*y(n,1))/0.11)+((20/0.11)-(6*(y(n+1,1)+y(n+1,2))+6*y(n+1,1))/0.11));
    y(n+1,2)=y(n,2)+(h/2)*(((20/0.3)-(6*(y(n,1)+y(n,2))+3*y(n,2))/0.3)+((20/0.3)-(6*(y(n+1,1)+y(n+1,2))+3*y(n+1,2))/0.3));  
    y(n+1,3)=y(n+1,1)+y(n+1,2);
end
plot(i1,y(:,3))
plot(i2,y(:,1),'r')
plot(i3,y(:,2),'g')
xlabel('Tiempo');
ylabel('Intensidad(t)');
hold off

Gráfica del resultado del Método de Euler y Trapecio en el cálculo de las intensidades.

En el resultado de la gráfica, lo obtenido con el método de Euler se representa con los símbolos "+", mientras que la línea continua corresponde al Método del Trapecio. Como se puede apreciar, ambos métodos dan resultados similares y muy próximos. En la interpretación de los resultados, diríamos que se trata de corrientes muy estables.

2.2 Variación de las intensidades de un circuito RL cuando se altera el valor de una de las resistencias.

Ahora vamos a analizar como afecta la variación del módulo de una resistencia (se incrementa, [math] R_3=3Ω [/math] pasa a valer [math] R_3=9Ω [/math]) en las intensidades del circuito. Nuestro análisis se realizará sustituyendo el nuevo valor en el mismo sistema, y se resolverá por el método de Euler y Trapecio comparando en una gráfica los resultados con los obtenidos anteriormente.

• MÉTODO DE EULER:

OCTAVE

%DATOS:
clc;clear all;
    %tiempo inicial
t0=0;
    %tiempo máximo a evaluar
tN=0.04;
    %soluciones en y0 para las variables i2 e i3
y0=0;
    %intervalos en los que se dividen
N=25;
	%salto de intervalo en intervalo
h=(tN-t0)/N;
    %vectores solución desde t0 a tN de h en h
i1=t0:h:tN;
i2=t0:h:tN;
i3=t0:h:tN;
    %vector de ceros y que va a almacenar las soluciones
    %tendra N+1 términos porque el primero es el cero y el ultimo el N
    %la intensidad dos es y(:,1)
    %la intensidad tres es y(:,2)
    %la intensidad uno es y(:,3)
y=zeros(N+1,3);
y(1,1)=y0;
y(1,2)=y0;
y(1,3)=y0;
%PARA RESOLVER POR EULER EL SISTEMA:
%y(n+1)=y(n)+h*f(t(n),y(n))
%f(t,i2)=((E/L2)-R1*(i2+i3)+R2*i2)/L2
%f(t,i3)=((E/L1)-R1*(i2+i3)+R3*i3)/L1
for n=1:N
    y(n+1,1)=y(n,1)+h*((20/0.11)-(6*(y(n,1)+y(n,2))+6*y(n,1))/0.11);
    y(n+1,2)=y(n,2)+h*((20/0.3)-(6*(y(n,1)+y(n,2))+9*y(n,2))/0.3);
    y(n+1,3)=y(n+1,1)+y(n+1,2);
end
hold on
plot(i1,y(:,3),'c+')
plot(i2,y(:,1),'m+')
plot(i3,y(:,2),'k+')

• MÉTODO DEL TRAPECIO:

OCTAVE

%DATOS:
clc;clear all;
    %tiempo inicial
t0=0;
    %tiempo máximo a evaluar
tN=0.04;
    %soluciones en y0 para las variables i2 e i3
y0=0;
    %intervalos en los que se dividen
N=25;
	%salto de intervalo en intervalo
h=(tN-t0)/N;
    %vectores solución desde t0 a tN de h en h
i1=t0:h:tN;
i2=t0:h:tN;
i3=t0:h:tN;
    %vector de ceros y que va a almacenar las soluciones
    %tendra N+1 términos porque el primero es el cero y el ultimo el N
    %la intensidad dos es y(:,1)
    %la intensidad tres es y(:,2)
    %la intensidad uno es y(:,3)
y=zeros(N+1,3);
y(1,1)=y0;
y(1,2)=y0;
y(1,3)=y0;
%PARA RESOLVER POR EL MÉTODO DEL TRAPECIO:
%y(n+1)=y(n)+(h/2)*(f(t(n),y(n)+f(t(n+1),y(n+1))
%f(t,i2)=((E/L2)-R1*(i2+i3)+R2*i2)/L2
%f(t,i3)=((E/L1)-R1*(i2+i3)+R3*i3)/L1
for n=1:N
    y(n+1,1)=y(n,1)+(h/2)*(((20/0.11)-(6*(y(n,1)+y(n,2))+6*y(n,1))/0.11)+((20/0.11)-(6*(y(n+1,1)+y(n+1,2))+6*y(n+1,1))/0.11));
    y(n+1,2)=y(n,2)+(h/2)*(((20/0.3)-(6*(y(n,1)+y(n,2))+9*y(n,2))/0.3)+((20/0.3)-(6*(y(n+1,1)+y(n+1,2))+9*y(n+1,2))/0.3));  
    y(n+1,3)=y(n+1,1)+y(n+1,2);
end
plot(i1,y(:,3),'c')
plot(i2,y(:,1),'m')
plot(i3,y(:,2),'k')
hold off

Gráfica de la intensidad del circuito RL compuesto tras incrementar una de sus resistencias

Gráfica del resultado del Método de Euler y Trapecio en el cálculo de las intensidades sin incrementar.

Gráfica comparativa de las intensidades de un circuito RL intensidades incrementadas y sin incrementar.

De nuevo, las cruces + corresponderían a lo obtenido por el método de Euler y la línea continua al Método del trapecio. En la gráfica superpuesta se observa el resultado final de la comparativa: frente a resistencias menores las intensidades resultan mayores, y viceversa. Esto se justifica gracias a la Ley de Ohm que establece [math] V=I*R[/math] de manera que la intensidad y resistencia son inversamente proporcionales, y si una se incrementa, la otra decrece: [math] I=V/R[/math], [math] R=V/I[/math].

3 Resolución de un problema de valor inicial en otro punto

3.1 Cálculo de las intensidades en un punto (el inicial) a partir de un valor inicial en un punto distinto al origen.

Esta vez vamos a calcular el valor concreto de un punto de una función, que hemos definido tras resolver una ecuación diferencial. Nuestra función será la que define las intensidades de un circuito RL compuesto. En la resolución de la ecuación diferencial que relaciona las derivadas con la función, procederemos de nuevo tratándolo como un problema de valor inicial; a pesar de que no sea un valor en el punto de origen el que tengamos por dato, sino un valor de [math] i_2 [/math] e [math] i_3 [/math] en el tiempo [math] t=0,02 [/math]. Serán precisamente estos valores en el origen los que calculemos. De nuevo, utilizaremos el segundo circuito expuesto en este artículo (sin incrementar resistencias):

Circuito RL compuesto

:[math]\frac{d}{dt}\cdot \begin{pmatrix} i_2\\i_3 \end{pmatrix} = - \begin{pmatrix} \frac{6+6}{0,11}&\frac{6}{0,11} \\ \frac{6}{0,3}&\frac{6+3}{0,3}\end{pmatrix} \cdot \begin{pmatrix} i_2\\i_3 \end{pmatrix} + \begin{pmatrix} \frac{20}{0,11} \\ \frac{20}{0,3} \end{pmatrix}[/math]

[math] R_1=R_2=6Ω [/math] , [math] R_3=3Ω [/math] , [math] L_1=0,3Hz [/math] , [math] L_2=0,11Hz [/math] y [math] E(t)=20 V [/math].

Nuestros valores iniciales [math] i_2(0,02)=i_3(0,02)=1A [/math]

Buscamos [math] i_1(t),i_2(t),i_3(t) [/math].

Resolvemos por el método de Euler:

%DATOS:
clc;clear all;
    %tiempo inicial
t0=0.02;
    %tiempo a evaluar
tN=0;
    %soluciones en y0.02 para las variables i2 e i3
y02=1;
    %
N=100000;
	%salto de intervalo en intervalo
h=(tN-t0)/N;
    %vectores solución desde t0 a tN de h en h
i1=t0:h:tN;
i2=t0:h:tN;
i3=t0:h:tN;
    %vector de ceros y que va a almacenar las soluciones
    %tendra N+1 términos porque el primero es el cero y el ultimo el N
    %la intensidad dos es y(:,1)
    %la intensidad tres es y(:,2)
    %la intensidad uno es y(:,3)
y=zeros(N+1,3);
y(1,1)=y02;
y(1,2)=y02;
y(1,3)=2;
%RESOLVEMOS POR EULER EL SISTEMA:
%y(n+1)=y(n)+h*f(t(n),y(n))
%f(t,i2)=((E/L2)-R1*(i2+i3)+R2*i2)/L2
%f(t,i3)=((E/L1)-R1*(i2+i3)+R3*i3)/L1
for n=1:N
    y(n+1,1)=y(n,1)+h*((20/0.11)-(6*(y(n,1)+y(n,2))+6*y(n,1))/0.11);
    y(n+1,2)=y(n,2)+h*((20/0.3)-(6*(y(n,1)+y(n,2))+3*y(n,2))/0.3);
    y(n+1,3)=y(n+1,1)+y(n+1,2);
end
intensidad1=y(N+1,3)
intensidad2=y(N+1,1)
intensidad3=y(N+1,2)

Los resultados obtenidos:

Valores de las intensidades para t=0

Así se prueba la eficacia de estas aproximaciones numéricas no sólo al trazar gráficas, sino también para dar valores concretos de la función en los puntos que se quiera.