MÉTODOS DE RUNGE – KUTTA

Introducción:
La computadora, es la herramienta mas poderosa hasta ahora conocida, para la solución de problemas en el campo de las ciencias exactas, en este caso los métodos numéricos, como punto principal por sus aplicaciones en la ingeniería.

Los métodos numéricos son técnicas, donde es posible resolver los problemas por medio de operaciones aritméticas, estos métodos implementan un buen numero de cálculos que son por demás demasiado lentos si se hacen manualmente, gastando mucha energía en la técnica misma de solución en vez de aplicarla sobre la definición del problema y su interpretación.

El trabajo monótono que se hacia anteriormente al uso de la computadora, hace de importancia, el dominio de los métodos numéricos, los cuales se deben llevar a cabo en combinación con las capacidades y potencialidades de la programación de computadoras para de esa forma resolver los problemas de ingeniería mucho mas fácilmente y eficientemente.


Objetivo de los métodos de Runge-Kutta:

El objetivo de los métodos numéricos de runge-kutta, es el análisis y solución de los problemas de valor inicial de ecuaciones diferenciales ordinarias (EDO), estos son una extensión del método de euler para resolver las (EDO’S), pero con un orden de exactitud mas alto que este.


Modelos matemáticos:

Los conocimientos científicos se usan rutinariamente por los ingenieros en el diseño de elementos tales como maquinas, circuitos eléctricos, estructuras etc.

Estos conocimientos son muy útiles cuando se expresan en forma de un modelo matemático, el cual se puede definir como una ecuación que expresa las características fundamentales de un sistema o proceso físico en términos matemáticos, siendo clasificados estos modelos, desde simples relaciones algebraicas hasta grandes y complicados sistemas de ecuaciones diferenciales.

Análisis de un modelo matemático:

Un modelo matemático: es una expresión matemática como veremos en el siguiente ejemplo:

Formula de la segunda ley de newton:

F= ma donde : F es la fuerza neta que actúa sobre el cuerpo.

m es la masa del objeto.

Utilizando esta ley, vamos a determinar la velocidad de un paracaidista en caída libre.

Para este caso puede crearse un nuevo modelo, expresando la aceleración como la razón de cambio de la velocidad con respecto al tiempo (dv/dt) .

Y sustituir en la ecuación de nueva forma:

F=m(dv/dt)

Así la masa, multiplicada por la razón de cambio de la velocidad es igual a la suma de fuerzas que actúan sobre el cuerpo.

Si la fuerza total cuando el objeto cae es positiva el objeto acelera, pero si es negativa desacelera, pero si la fuerza neta es cero la velocidad permanecerá a un nivel constante.

Para un cuerpo que cae la fuerza , la fuerza total esta compuesta por dos fuerzas contrarias, la atracción debida a la gravedad Fd, y la fuerza hacia arriba debida a la resistencia del aire Fu.

Por lo tanto: F= Fd + Fu

La fuerza debida a la gravedad Fd se puede reescribir:

Fd=mg donde g es la constante de gravitación que equivale a 980 cm por segundo cuadrado.

La resistencia del aire se puede formular como una aproximación sencilla linealmente proporcional a la velocidad:

Fv = -cv donde c es una constante de proporcionalidad llamada coeficiente de arrastre.

Entonces la fuerza total es la diferencia de las fuerzas hacia abajo y hacia arriba, así que combinando las ecuaciones anteriores:

M(dv/dt)= mg – cv o dividiendo cada lado entre m:

esta es la ecuación de un cuerpo que cae a las fuerzas que actúan sobre el y es una ecuación diferencial porque esta escrita en términos de la razón diferencial dv/dt.

Usando el calculo y resolviendo esta ecuación diferencial se puede llegar a la siguiente formula que expresa la velocidad del paracaidista en función del tiempo.


 Ec(principal).

A continuación se observara la diferencia existente entre tres tipos de soluciones para el problema del paracaidista analizando las ventajas de cada uno de ellos con respecto a los demás:

Enunciado del problema de un paracaidista que cae:


MÉTODO ANALÍTICO:

Un paracaidista con una masa de 55500 g salta de un aeroplano apliquese la ecuación principal para calcular la velocidad antes de abrir el paracaídas. El cociente de arrastre c es aproximadamente igual a 10500 g/s.

Solución :

Al sustituir los valores de los parámetros en la ecuación principal se obtiene

Al dar varios valores de t se obtienen las velocidades se obtienen las velocidades para el tiempo, los resultados se presentan a continuación:

TABLA DE RESULTADOS

Tiempo en segundos Velocidad en cm/s.
0 0
2 1631.7
4 2749.5
6 3515.1
8 4039.6
10 4398.8
12 4644.9
Al infinito 5180.0

Nota: La escala de la velocidad en la gráfica es de 1=1000.

Si lo has notado para sacar el resultado de tus cálculos de la tabla y gráfica anteriores necesitas estar sustituyendo en la formula de v(t) esto hace el método analítico cansado y repetitivo, pero es una solución analítica exacta porque satisface la ecuación diferencial original.


SOLUCIÓN NUMÉRICA:

La aproximación de la razón de cambio de la velocidad con respecto al tiempo puede representarse de la siguiente forma

Y llegamos a ordenar esta ecuación de la siguiente manera

la cual puede emplearse para extender el calculo , así que aplicandola directamente se dirá

Nuevo valor = valor anterior + valor estimulado x incremento del tiempo.

de v. de v de la pendiente

Así pues efectuando el mismo calculo para el problema se procede como sigue:
Sustituyendo en la ecuación, nuevamente vemos que necesitamos aplicar las matemáticas sin el auxilio de una computadora, y nos vamos a tardar un buen rato .

v(4)= 2330.8 cm/s.

Y así sucesivamente se continua con el calculo para obtener los valores en la siguiente tabla
 

Tiempo en seg. Vel. En cm/s.
0  
2  
4  
6  
8  
10  
12  
al infinito.   

Puede verse por el resultado en la tabla anterior que la solución por un método numérico se aproxima bastante bien a la solución exacta, pero debido al empleo de rectas para aproximar la función que es continuamente curva, existe discrepancia entre los resultados de la tabla del método analítico y la de este método , una manera de minimizar el error es utilizando intervalos menores de tiempo en el muestreo de la ecuación, y así los segmentos de recta seguirán mas de cerca a la solución verdadera.

Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequeños hará poco practicas las soluciones numéricas, es por eso que a continuación se presentara un método con el cual se hará uso de un programa de computadora para la solución del problema.
 


SOLUCIÓN POR EL MÉTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA.

Los métodos de euler son procedimientos para resolver EDO’S de primer orden y se pueden programar en una computadora para no hacer los cálculos manualmente.

Una desventaja importante de los métodos de euler es que el orden de exactitud es bajo comparado con los de RUNGE- KUTTA .Y aunque para uno y otro existen versiones diferentes nosotros solamente nos basaremos en su potencialidad de aplicación dentro de la programación por computadora para la solución del problemas.


PROGRAMA PARA LA SOLUCIÓN DEL PROBLEMA DEL PARACAIDISTA.

clear, clf, hold off
t=0; n=0; v=0;
C=10.5; M=55.5; g=9.8; h=.1;
t_rec(1)=t; v_rec(1)=v;
while t<=12
n=n+1;
v=v+h*(-C/M*v*v+g);
t=t+h;
v_rec(n+1)=v;
t_rec(n+1)=t;
end
plot(t_rec,v_rec)

Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programación como fortran o C, pero mas sencillo que estos puesto que cualquier variable puede contener números de cualquier tipo sin una declaración especial durante la programación .

Por lo que resulta una poderosa herramienta para el análisis numérico y matemático en la aplicación de los métodos numéricos a la solución de los problemas de ingeniería.
La siguiente gráfica muestra el resultado del programa anterior (también realizada por matlab).

Velocidad cm/s.
Tiempo en seg.

Podemos observar, con los resultados anteriores , que la solución por medio de la programación por computadora se puede hacer que se aproxime a la solución exacta de la ecuación diferencial tanto como se quiera, dependiendo del método que se utilice.

En los métodos de RUNGE-KUTTA, el método de exactitud se incrementa mediante el empleo de un método de integración numérica de mas alto orden, la mayor exactitud implica que el resultado calculado es mas exacto y que los errores se reducen con mayor rapidez al reducirse h, donde h es un intervalo de tiempo fijo que se utiliza repetidamente.

Si consideramos la ecuación diferencial ordinaria siguiente

Si queremos calcular

Aplicamos una regla trapezoidal al miembro derecho , la cual resulta de dividir el área de una curva en n intervalos iguales de longitud h evaluando el área de cada una de las fajas, aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita, si h es lo suficientemente pequeña, entonces se puede aproximar como el área de un trapecio bajo la curva, es decir:

El área total es la suma de las áreas parciales

Entonces la forma trapezoidal en donde el área se aproxima mediante una suma de n trapecios es

Esta formula es muy importante ya que se aplica directamente a las diferentes versiones de los métodos de RUNGE-KUTTA , para poder aplicar las ecuaciones a problemas concretos.

A continuación veremos algunos ejemplos de los métodos de runge a problemas de circuitos eléctricos, y su solución mas rápida y eficiente por medio de un programa realizado en computadora por MATLAB que como ya vimos anteriormente , ahorra la fatiga de calcular manualmente.

El comportamiento de un circuito eléctrico cambia significativamente dependiendo de los valores de los componentes empleados , así, en el circuito que se muestra a continuación , la inductancia L=50mH, una resistencia R=20ohms y una fuente de voltaje de E=10v. entonces, si se cierra el interruptor en un t=0 la corriente I(t) satisface la ecuación diferencial

I(0) = 0

Se necesita encontrar el valor de la corriente para 0< t < 0.02 s. mediante un método de RUNGE-KUTTA de segundo orden con h=.0001.


donde

Mostramos entonces los cálculos para los dos primeros pasos:

t = 0

t= 0.0001

en el segundo paso se procede de la misma manera, obteniendo para la corriente dos un resultado de .038431, pero como podemos darnos cuenta , seria muy tardado resolverlo manualmente de ahí que lo que falta para calcular se hará en el programa de matlab.

NOTA: Los programas pueden estar sujetos a errores al escribirse, así que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas fácilmente.


PROGRAMA PARA RESOLVER EL CIRCUITO ANTERIOR:
clear, clf, hold off
R=20; %ohm
L=50e-3; %H
E=10; %V
y(1)=0; t(1)=0;
h=.1e-3;
n=1;
y_rec(1)=y; t_rec(1)=0; t=0;
RL=R/L; EL=E/L;
while t(n)<.02
K1=h*fn10_9(y(n),RL,EL);
K2=h*fn10_9(y(n)+K1,RL,EL);
y(n+1)=y(n)+.5*(K1+K2);
t(n+1)=n+h;
n=n+1;
end
plot(t,y)
xlabel('t = seg. ' )
ylabel('I(A)')
fn10_9
function f=fn10_9(I,RL,EL)
f=-RL*I+EL

Ahora vamos a un segundo ejemplo de aplicación a circuitos y veamos la siguiente figura

Los valores de los elementos que se encuentran en la figura anterior son hipotéticos es solo para representar el circuito, el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solución por medio de un listado de programa.
Para la primera inductancia, resistencia y capacitor

con
donde e(t) = 0 excepto e(t)=1 cuando 0<t<0.01s, q(t) es la carga del condensador, i1(t) e i2(t) son corrientes; Las condiciones iniciales son i1(0) = i2(0) = q(0) = 0 en todos los casos.

Resolviendo las ecuaciones:

 a) 
b)igual que (a) excepto que
 c)igual que (a) excepto que
 d)igual que (a) excepto que
La solución se determina en el siguiente programa realizado en matlab

clear;clg
subplot(221)
for k=1:4
e=1;
if k=1; subplot (221);
La=.01; Lb=.5; Ra=200; Rb=20; C=.002; end
if k=2 subplot(222);
La=.1; Lb=.5; Ra=200; Rb=20;C=.002; end
if k=3; subplot(223);
La=.01; Lb=.25; Ra=200; Rb=20; C=.002; end
if k=4; subplot(224)
La=.01; Lb=.5; Ra=20; Rb=20; C=.002; end
M=[-Ra/La, Ra/La, -1/(La*C);...
Ra/Lb,-(Ra+Rb)/Lb, 1/(Lb*C);...
1/C,-1/C,0]
S=[0;0;0]; X=[0;0;0];
h=.00005;
for n=1:101
t=(n-1)*h;
%S=[sin(t*600)*exp(-t*600)/La;0;0]
%S=[cos(t*600)/La;0;0];
S=[1/La;0;0];
if t>.001, S=[0;0;0];end
k1=h*(M*X+S);
k2=h*(M*(X+k1/2)+S);
k3=h*(M*(X+k2/2)+S);
k4=h*(M*(X+k3)+S);
X=X+(K1+K2*2+k3*2+k4)/6;
X_r(:,n)=X;
t_r(n)=t;
end

plot(t_r,X_r(1:2,:),t_r,X_r(3,:))
xlabel('t'),ylabel('i1,i2,q')
L=length(t_r)
text(t_r(L/10),X_r(1,L/10),'i1')
text(t_r(L/2),X_r(2,L/2),'i2')
text(t_r(L*.8),X_r(3,L*.8),'q')
if k==1;title('(A)'),end
if k==2;title('(B)');end
if k==3;title('(C)');end
if k==4;title('(D)'),end
end

Hay que tener en cuenta que de alguna u otra manera los métodos numéricos estan sujetos a ciertos errores al hacer los cálculos en muchos pasos (iteraciones), pero que estos se pueden llegar a despreciar por la actual eficiencia de las maquinas computadoras, las cuales pueden llegar a tomar los intervalos de tiempo infinitesimalmente pequeños para la aplicación de las ecuaciones y obtener los resultados mas rápido y de manera mas exacta.

No hay limite para las aplicaciones que se pueden realizar en programación de computadoras para resolver los problemas de métodos numéricos así que ya para terminar a continuación se darán algunos otros ejemplos en matlab para algunos métodos mencionados como clásicos en algunos libros .


Programa uno: function.
function y = f1_(x)
y = 0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5;


Programa dos: Que toma al programa 1 como archivo para funcionar (una de las características mas importantes de matlab.).

Método de Newton:
function newton_(x0)


%SOLUCION DE ECUACIONES NO LINEALES

%27 DE NOVIEMBRE DE 1997
%"METODO DE ITERACION DE NEWTON"
% LA FUNCION ESTA EN "f1_.m"
% "x0" estimación incial
x=x0; xb=x-999;
n=0; dx=0.01;
while abs(x-xb)>0.000001
n=n+1; xb=x;
if n>300 break; end
y=f1_(x);
ydr=(f1_(x+dx)-y)/dx;
x=xb-y/ydr;
fprintf('n=%3.0f, x=%12.5e, y=%12.5e, ', n,x,y)
fprintf('yd = %12.5e \n', ydr)
end
fprintf('\n Respuesta final = %12.6e\n',x);


Programa tres: Metodo de la regla falsa.

function rfalsa_(a,c)
SOLUCIÓN DE ECUACIONES NO LINEALES
%23 DE NOVIEMBRE DE 1997
%"METODO DE REGLA FALSA"
% LA FUNCION ESTA EN "f1_.m"
% "a" y "c" intervalo donde se cree esta la raiz
tolerancia=0.000001; iteraciones=30;
fprintf('\n it a b c f(a) f(b) f(c) \n');
it=0;
ya=f1_(a); yc=f1_(c);
if (ya*yc > 0)
fprintf('\n\n no hay raiz en este intervalo \n');
else
while 1
it=it+1;
b=(a+c)/2; yb=f1_(b);
fprintf('%3.0f %6.5f, %6.5f %10.5f %10.5f %10.5f \n',it,a,b,c,ya,yb,yc);
if (abs(c-a)/2 <= tolerancia)
fprintf('se cumplio la tolerancia. \n'); break
fprintf('\n Cambie de intervalo y ejecute de nuevo. \n');
end
if (it>iteraciones)
fprintf('se exede el numero de iteraciones.\n');break
end
if (ya*yb<=0) c=b; yc=yb;
else a=b; ya=yb;
end
end
fprintf('\n Raíz=%12.6f\n',b);
end


Programa cuatro: Metodo de simpson. 1/3.

%METODOS DE INTEGRACION
%18 DE NOVIEMBRE DE 1997
%"METODO DE SIMPSON 1/3"
% LA FUNCION A INTEGAR ESTA EN "f1_.m"
Iexacta=1.64053334;
a=0; b=0.8; n=20;
h=(b-a)/n; c=1:n+1;
x=a+(c-1)*h;
Is=(h/3)*(f1_(x(1))+4*sum(f1_(x(2:2:n)))+f1_(x(n+1)));
if n>2 Is=Is+(h/3)*2*sum(f1_(x(3:2:n))); end
error=(Is/Iexacta)*100;
fprintf('\n método de simpson 1/3 \n');
fprintf('\n n Is Iexacta exactitud\n');
fprintf('%3.0f %10.5f %10.5f %10.5f\n', n, Is, Iexacta, error);


Programa cinco: Metodo de simpson 3/8.

%METODOS DE INTEGRACION
%18 DE NOVIEMBRE DE 1997
%"METODO DE SIMPSON 3/8"
% LA FUNCION A INTEGAR ESTA EN "f1_.m"
Iexacta=1.64053334;
a=0; b=0.8; n=3;
h=(b-a)/n; c=1:n+1;
x=a+(c-1)*h;
Is=(3/8)*h*(f1_(x(n-2))+3*(f1_(x(n-1)))+3*(f1_(x(n)))+f1_(x(n+1)));
if n>3 Is=Is+(h/3)*2*sum(f1_(x(3:2:n))); end
error=(Is/Iexacta)*100;
fprintf('\n método de simpson 1/3 \n');
fprintf('\n n Is Iexacta exactitud\n');
fprintf('%3.0f %10.5f %10.5f %10.5f\n', n, Is, Iexacta, error);
 


Programa seis: Metodo de la regla trapezoidal .

%METODOS DE INTEGRACION
%18 DE NOVIEMBRE DE 1997
%"METODO TRAPEZOIDAL"
% LA FUNCION A INTEGAR ESTA EN "f1_.m"
Iexacta=1.64053334;
a=0; b=0.8; n=20;
h=(b-a)/n;
x=a:h:b;
It=(h/2)*(2*sum(f1_(x(2:n)))+f1_(x(1))+f1_(x(n+1)));
exactitud=(It/Iexacta)*100;
fprintf('\n método trapezoidal \n');
fprintf('\n n It Iexacta exactitud\n');
fprintf('%3.0f %10.5f %10.5f %10.5f\n', n, It, Iexacta, exactitud);


BIBLIOGRAFIA:

-Apuntes de la maestría en ciencias de ingenieria electronica, UdeG..
-Calculo diferencial e integral. William Anthony Granville. Ed. LIMUSA.
-Calculo con Geometria Analitica. Earl W. Swokowsky. Grupo Editorial Iberoamerica.
-Instrumentacion electronica. a. j. Diefenderfer. Interamericana.
-El lenguaje de programación C. Brian W. Kernighan. Prentice Hall.
-Análisis numérico y visualizacion gráfica con matlab.. Soichiro Nakamura. Prentice Hall.
-Metodos Numericos para Ingenieros con aplicaciones a computadoras personales.
Steven C Chapra. Prentice Hall.