Post on 10-Jun-2018
transcript
UNIDAD 5- Solucioacuten de ecuaciones diferenciales
51 Meacutetodos de un paso
511 Meacutetodo de Euler y Euler mejorado
Meacutetodo de Euler
Se llama meacutetodo de Euler al meacutetodo numeacuterico consistente en ir incrementando paso a paso la variable independiente y hallando la siguiente imagen con la derivada
La primera derivada proporciona una estimacioacuten directa de la pendiente en Xi (ver Graacutefico
Nordm01) [1]
Donde f (Xi Yi) es la ecuacioacuten diferencial evaluada en Xi y Yi Tal estimacioacuten podraacute substituirse en la ecuacioacuten [2] nos queda que
[2]
Esta foacutermula es conocida como el meacutetodo de Euler (punto medio) Se predice un nuevo valor de Y por medio de la pendiente (igual a la primera derivada en el valor original de X)
Error para el meacutetodo de Euler
La solucioacuten numeacuterica de las ecuaciones diferenciales ordinarias (EDO) involucra dos tipos de error
1) Errores de Truncamiento o discretizacion causados por la naturaleza de las teacutecnicas empleadas para aproximar los valores de y
2) Errores de Redondeo que son el resultado del nuacutemero limite de cifras significativas que pueden retener una computadora
Meacutetodo de Euler Mejorado
Este meacutetodo se basa en la misma idea del meacutetodo anterior pero hace un refinamiento en la aproximacioacuten tomando un promedio entre ciertas pendientes
La foacutermula es la siguiente
Donde
Para entender esta foacutermula analicemos el primer paso de la aproximacioacuten con base en la siguiente graacutefica
En la graacutefica vemos que la pendiente promedio corresponde a la pendiente de la recta bisectriz de la recta tangente a la curva en el punto de la condicioacuten inicial y la recta
tangente a la curva en el punto donde es la aproximacioacuten obtenida con la primera foacutermula de Euler Finalmente esta recta bisectriz se traslada paralelamente hasta
el punto de la condicioacuten inicial y se considera el valor de esta recta en el punto como la aproximacioacuten de Euler mejorada
512 Meacutetodo de Runge-Kutta
Introduccioacuten La computadora es la herramienta mas poderosa hasta ahora conocida para la solucioacuten de problemas en el campo de las ciencias exactas en este caso los meacutetodos numeacutericos como punto principal por sus aplicaciones en la ingenieriacutea
Los meacutetodos numeacutericos son teacutecnicas donde es posible resolver los problemas por medio de operaciones aritmeacuteticas estos meacutetodos implementan un buen numero de caacutelculos que son por demaacutes demasiado lentos si se hacen manualmente gastando mucha energiacutea en la teacutecnica misma de solucioacuten en vez de aplicarla sobre la definicioacuten del problema y su interpretacioacuten
El trabajo monoacutetono que se hacia anteriormente al uso de la computadora hace de importancia el dominio de los meacutetodos numeacutericos los cuales se deben llevar a cabo en combinacioacuten con las capacidades y potencialidades de la programacioacuten de computadoras para de esa forma resolver los problemas de ingenieriacutea mucho mas faacutecilmente y eficientemente
Objetivo de los meacutetodos de Runge-Kutta
El objetivo de los meacutetodos numeacutericos de runge-kutta es el anaacutelisis y solucioacuten de los problemas de valor inicial de ecuaciones diferenciales ordinarias (EDO) estos son una extensioacuten del meacutetodo de euler para resolver las (EDOrsquoS) pero con un orden de exactitud mas alto que este
Modelos matemaacuteticos
Los conocimientos cientiacuteficos se usan rutinariamente por los ingenieros en el disentildeo de elementos tales como maquinas circuitos eleacutectricos estructuras etc
Estos conocimientos son muy uacutetiles cuando se expresan en forma de un modelo matemaacutetico el cual se puede definir como una ecuacioacuten que expresa las caracteriacutesticas fundamentales de un sistema o proceso fiacutesico en teacuterminos matemaacuteticos siendo clasificados estos modelos desde simples relaciones algebraicas hasta grandes y complicados sistemas de ecuaciones diferenciales
Anaacutelisis de un modelo matemaacutetico
Un modelo matemaacutetico es una expresioacuten matemaacutetica como veremos en el siguiente ejemplo
Formula de la segunda ley de newton
F= ma donde F es la fuerza neta que actuacutea sobre el cuerpo
m es la masa del objeto
Utilizando esta ley vamos a determinar la velocidad de un paracaidista en caiacuteda libre
Para este caso puede crearse un nuevo modelo expresando la aceleracioacuten como la razoacuten de cambio de la velocidad con respecto al tiempo (dvdt)
Y sustituir en la ecuacioacuten de nueva forma
F=m(dvdt)
Asiacute la masa multiplicada por la razoacuten de cambio de la velocidad es igual a la suma de fuerzas que actuacutean 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 permaneceraacute a un nivel constante
Para un cuerpo que cae la fuerza la fuerza total esta compuesta por dos fuerzas contrarias la atraccioacuten 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 gravitacioacuten que equivale a 980 cm por segundo cuadrado
La resistencia del aire se puede formular como una aproximacioacuten 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 asiacute que combinando las ecuaciones anteriores
M(dvdt)= mg ndash cv o dividiendo cada lado entre m
esta es la ecuacioacuten de un cuerpo que cae a las fuerzas que actuacutean sobre el y es una ecuacioacuten diferencial porque esta escrita en teacuterminos de la razoacuten diferencial dvdt
Usando el calculo y resolviendo esta ecuacioacuten diferencial se puede llegar a la siguiente formula que expresa la velocidad del paracaidista en funcioacuten del tiempo
Ec(principal)
A continuacioacuten 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 demaacutes
Enunciado del problema de un paracaidista que cae
MEacuteTODO ANALIacuteTICO
Un paracaidista con una masa de 55500 g salta de un aeroplano apliquese la ecuacioacuten principal para calcular la velocidad antes de abrir el paracaiacutedas El cociente de arrastre c es aproximadamente igual a 10500 gs
Solucioacuten
Al sustituir los valores de los paraacutemetros en la ecuacioacuten 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 continuacioacuten
TABLA DE RESULTADOS
Tiempo en segundos
Velocidad en cms
0 0
2 16317
4 27495
6 35151
8 40396
10 43988
12 46449
Al infinito 51800
Nota La escala de la velocidad en la graacutefica es de 1=1000
Si lo has notado para sacar el resultado de tus caacutelculos de la tabla y graacutefica anteriores necesitas estar sustituyendo en la formula de v(t) esto hace el meacutetodo analiacutetico cansado y repetitivo pero es una solucioacuten analiacutetica exacta porque satisface la ecuacioacuten diferencial original
SOLUCIOacuteN NUMEacuteRICA
La aproximacioacuten de la razoacuten de cambio de la velocidad con respecto al tiempo puede representarse de la siguiente forma
Y llegamos a ordenar esta ecuacioacuten de la siguiente manera
la cual puede emplearse para extender el calculo asiacute que aplicandola directamente se diraacute
Nuevo valor = valor anterior + valor estimulado x incremento del tiempo
de v de v de la pendiente
Asiacute pues efectuando el mismo calculo para el problema se procede como sigue Sustituyendo en la ecuacioacuten nuevamente vemos que necesitamos aplicar las matemaacuteticas sin el auxilio de una computadora y nos vamos a tardar un buen rato
v(4)= 23308 cms
Y asiacute sucesivamente se continua con el calculo para obtener los valores en la siguiente tabla
Tiempo en seg
Vel En cms
0
2
4
6
8
10
12
al infinito
Puede verse por el resultado en la tabla anterior que la solucioacuten por un meacutetodo numeacuterico se aproxima bastante bien a la solucioacuten exacta pero debido al empleo de rectas para aproximar la funcioacuten que es continuamente curva existe discrepancia entre los resultados de la tabla del meacutetodo analiacutetico y la de este meacutetodo una manera de minimizar el error es utilizando intervalos menores de
tiempo en el muestreo de la ecuacioacuten y asiacute los segmentos de recta seguiraacuten mas de cerca a la solucioacuten verdadera
Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequentildeos haraacute poco practicas las soluciones numeacutericas es por eso que a continuacioacuten se presentara un meacutetodo con el cual se haraacute uso de un programa de computadora para la solucioacuten del problema
SOLUCIOacuteN POR EL MEacuteTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA
Los meacutetodos de euler son procedimientos para resolver EDOrsquoS de primer orden y se pueden programar en una computadora para no hacer los caacutelculos manualmente
Una desventaja importante de los meacutetodos 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 aplicacioacuten dentro de la programacioacuten por computadora para la solucioacuten del problemas
PROGRAMA PARA LA SOLUCIOacuteN DEL PROBLEMA DEL PARACAIDISTA
clear clf hold off t=0 n=0 v=0 C=105 M=555 g=98 h=1 t_rec(1)=t v_rec(1)=v while tlt=12 n=n+1 v=v+h(-CMvv+g) t=t+h v_rec(n+1)=v t_rec(n+1)=t end plot(t_recv_rec)
Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programacioacuten como fortran o C pero mas sencillo que estos puesto que cualquier variable puede contener nuacutemeros de cualquier tipo sin una declaracioacuten especial durante la programacioacuten
Por lo que resulta una poderosa herramienta para el anaacutelisis numeacuterico y matemaacutetico en la aplicacioacuten de los meacutetodos numeacutericos a la solucioacuten de los problemas de ingenieriacutea
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
511 Meacutetodo de Euler y Euler mejorado
Meacutetodo de Euler
Se llama meacutetodo de Euler al meacutetodo numeacuterico consistente en ir incrementando paso a paso la variable independiente y hallando la siguiente imagen con la derivada
La primera derivada proporciona una estimacioacuten directa de la pendiente en Xi (ver Graacutefico
Nordm01) [1]
Donde f (Xi Yi) es la ecuacioacuten diferencial evaluada en Xi y Yi Tal estimacioacuten podraacute substituirse en la ecuacioacuten [2] nos queda que
[2]
Esta foacutermula es conocida como el meacutetodo de Euler (punto medio) Se predice un nuevo valor de Y por medio de la pendiente (igual a la primera derivada en el valor original de X)
Error para el meacutetodo de Euler
La solucioacuten numeacuterica de las ecuaciones diferenciales ordinarias (EDO) involucra dos tipos de error
1) Errores de Truncamiento o discretizacion causados por la naturaleza de las teacutecnicas empleadas para aproximar los valores de y
2) Errores de Redondeo que son el resultado del nuacutemero limite de cifras significativas que pueden retener una computadora
Meacutetodo de Euler Mejorado
Este meacutetodo se basa en la misma idea del meacutetodo anterior pero hace un refinamiento en la aproximacioacuten tomando un promedio entre ciertas pendientes
La foacutermula es la siguiente
Donde
Para entender esta foacutermula analicemos el primer paso de la aproximacioacuten con base en la siguiente graacutefica
En la graacutefica vemos que la pendiente promedio corresponde a la pendiente de la recta bisectriz de la recta tangente a la curva en el punto de la condicioacuten inicial y la recta
tangente a la curva en el punto donde es la aproximacioacuten obtenida con la primera foacutermula de Euler Finalmente esta recta bisectriz se traslada paralelamente hasta
el punto de la condicioacuten inicial y se considera el valor de esta recta en el punto como la aproximacioacuten de Euler mejorada
512 Meacutetodo de Runge-Kutta
Introduccioacuten La computadora es la herramienta mas poderosa hasta ahora conocida para la solucioacuten de problemas en el campo de las ciencias exactas en este caso los meacutetodos numeacutericos como punto principal por sus aplicaciones en la ingenieriacutea
Los meacutetodos numeacutericos son teacutecnicas donde es posible resolver los problemas por medio de operaciones aritmeacuteticas estos meacutetodos implementan un buen numero de caacutelculos que son por demaacutes demasiado lentos si se hacen manualmente gastando mucha energiacutea en la teacutecnica misma de solucioacuten en vez de aplicarla sobre la definicioacuten del problema y su interpretacioacuten
El trabajo monoacutetono que se hacia anteriormente al uso de la computadora hace de importancia el dominio de los meacutetodos numeacutericos los cuales se deben llevar a cabo en combinacioacuten con las capacidades y potencialidades de la programacioacuten de computadoras para de esa forma resolver los problemas de ingenieriacutea mucho mas faacutecilmente y eficientemente
Objetivo de los meacutetodos de Runge-Kutta
El objetivo de los meacutetodos numeacutericos de runge-kutta es el anaacutelisis y solucioacuten de los problemas de valor inicial de ecuaciones diferenciales ordinarias (EDO) estos son una extensioacuten del meacutetodo de euler para resolver las (EDOrsquoS) pero con un orden de exactitud mas alto que este
Modelos matemaacuteticos
Los conocimientos cientiacuteficos se usan rutinariamente por los ingenieros en el disentildeo de elementos tales como maquinas circuitos eleacutectricos estructuras etc
Estos conocimientos son muy uacutetiles cuando se expresan en forma de un modelo matemaacutetico el cual se puede definir como una ecuacioacuten que expresa las caracteriacutesticas fundamentales de un sistema o proceso fiacutesico en teacuterminos matemaacuteticos siendo clasificados estos modelos desde simples relaciones algebraicas hasta grandes y complicados sistemas de ecuaciones diferenciales
Anaacutelisis de un modelo matemaacutetico
Un modelo matemaacutetico es una expresioacuten matemaacutetica como veremos en el siguiente ejemplo
Formula de la segunda ley de newton
F= ma donde F es la fuerza neta que actuacutea sobre el cuerpo
m es la masa del objeto
Utilizando esta ley vamos a determinar la velocidad de un paracaidista en caiacuteda libre
Para este caso puede crearse un nuevo modelo expresando la aceleracioacuten como la razoacuten de cambio de la velocidad con respecto al tiempo (dvdt)
Y sustituir en la ecuacioacuten de nueva forma
F=m(dvdt)
Asiacute la masa multiplicada por la razoacuten de cambio de la velocidad es igual a la suma de fuerzas que actuacutean 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 permaneceraacute a un nivel constante
Para un cuerpo que cae la fuerza la fuerza total esta compuesta por dos fuerzas contrarias la atraccioacuten 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 gravitacioacuten que equivale a 980 cm por segundo cuadrado
La resistencia del aire se puede formular como una aproximacioacuten 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 asiacute que combinando las ecuaciones anteriores
M(dvdt)= mg ndash cv o dividiendo cada lado entre m
esta es la ecuacioacuten de un cuerpo que cae a las fuerzas que actuacutean sobre el y es una ecuacioacuten diferencial porque esta escrita en teacuterminos de la razoacuten diferencial dvdt
Usando el calculo y resolviendo esta ecuacioacuten diferencial se puede llegar a la siguiente formula que expresa la velocidad del paracaidista en funcioacuten del tiempo
Ec(principal)
A continuacioacuten 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 demaacutes
Enunciado del problema de un paracaidista que cae
MEacuteTODO ANALIacuteTICO
Un paracaidista con una masa de 55500 g salta de un aeroplano apliquese la ecuacioacuten principal para calcular la velocidad antes de abrir el paracaiacutedas El cociente de arrastre c es aproximadamente igual a 10500 gs
Solucioacuten
Al sustituir los valores de los paraacutemetros en la ecuacioacuten 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 continuacioacuten
TABLA DE RESULTADOS
Tiempo en segundos
Velocidad en cms
0 0
2 16317
4 27495
6 35151
8 40396
10 43988
12 46449
Al infinito 51800
Nota La escala de la velocidad en la graacutefica es de 1=1000
Si lo has notado para sacar el resultado de tus caacutelculos de la tabla y graacutefica anteriores necesitas estar sustituyendo en la formula de v(t) esto hace el meacutetodo analiacutetico cansado y repetitivo pero es una solucioacuten analiacutetica exacta porque satisface la ecuacioacuten diferencial original
SOLUCIOacuteN NUMEacuteRICA
La aproximacioacuten de la razoacuten de cambio de la velocidad con respecto al tiempo puede representarse de la siguiente forma
Y llegamos a ordenar esta ecuacioacuten de la siguiente manera
la cual puede emplearse para extender el calculo asiacute que aplicandola directamente se diraacute
Nuevo valor = valor anterior + valor estimulado x incremento del tiempo
de v de v de la pendiente
Asiacute pues efectuando el mismo calculo para el problema se procede como sigue Sustituyendo en la ecuacioacuten nuevamente vemos que necesitamos aplicar las matemaacuteticas sin el auxilio de una computadora y nos vamos a tardar un buen rato
v(4)= 23308 cms
Y asiacute sucesivamente se continua con el calculo para obtener los valores en la siguiente tabla
Tiempo en seg
Vel En cms
0
2
4
6
8
10
12
al infinito
Puede verse por el resultado en la tabla anterior que la solucioacuten por un meacutetodo numeacuterico se aproxima bastante bien a la solucioacuten exacta pero debido al empleo de rectas para aproximar la funcioacuten que es continuamente curva existe discrepancia entre los resultados de la tabla del meacutetodo analiacutetico y la de este meacutetodo una manera de minimizar el error es utilizando intervalos menores de
tiempo en el muestreo de la ecuacioacuten y asiacute los segmentos de recta seguiraacuten mas de cerca a la solucioacuten verdadera
Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequentildeos haraacute poco practicas las soluciones numeacutericas es por eso que a continuacioacuten se presentara un meacutetodo con el cual se haraacute uso de un programa de computadora para la solucioacuten del problema
SOLUCIOacuteN POR EL MEacuteTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA
Los meacutetodos de euler son procedimientos para resolver EDOrsquoS de primer orden y se pueden programar en una computadora para no hacer los caacutelculos manualmente
Una desventaja importante de los meacutetodos 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 aplicacioacuten dentro de la programacioacuten por computadora para la solucioacuten del problemas
PROGRAMA PARA LA SOLUCIOacuteN DEL PROBLEMA DEL PARACAIDISTA
clear clf hold off t=0 n=0 v=0 C=105 M=555 g=98 h=1 t_rec(1)=t v_rec(1)=v while tlt=12 n=n+1 v=v+h(-CMvv+g) t=t+h v_rec(n+1)=v t_rec(n+1)=t end plot(t_recv_rec)
Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programacioacuten como fortran o C pero mas sencillo que estos puesto que cualquier variable puede contener nuacutemeros de cualquier tipo sin una declaracioacuten especial durante la programacioacuten
Por lo que resulta una poderosa herramienta para el anaacutelisis numeacuterico y matemaacutetico en la aplicacioacuten de los meacutetodos numeacutericos a la solucioacuten de los problemas de ingenieriacutea
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
2) Errores de Redondeo que son el resultado del nuacutemero limite de cifras significativas que pueden retener una computadora
Meacutetodo de Euler Mejorado
Este meacutetodo se basa en la misma idea del meacutetodo anterior pero hace un refinamiento en la aproximacioacuten tomando un promedio entre ciertas pendientes
La foacutermula es la siguiente
Donde
Para entender esta foacutermula analicemos el primer paso de la aproximacioacuten con base en la siguiente graacutefica
En la graacutefica vemos que la pendiente promedio corresponde a la pendiente de la recta bisectriz de la recta tangente a la curva en el punto de la condicioacuten inicial y la recta
tangente a la curva en el punto donde es la aproximacioacuten obtenida con la primera foacutermula de Euler Finalmente esta recta bisectriz se traslada paralelamente hasta
el punto de la condicioacuten inicial y se considera el valor de esta recta en el punto como la aproximacioacuten de Euler mejorada
512 Meacutetodo de Runge-Kutta
Introduccioacuten La computadora es la herramienta mas poderosa hasta ahora conocida para la solucioacuten de problemas en el campo de las ciencias exactas en este caso los meacutetodos numeacutericos como punto principal por sus aplicaciones en la ingenieriacutea
Los meacutetodos numeacutericos son teacutecnicas donde es posible resolver los problemas por medio de operaciones aritmeacuteticas estos meacutetodos implementan un buen numero de caacutelculos que son por demaacutes demasiado lentos si se hacen manualmente gastando mucha energiacutea en la teacutecnica misma de solucioacuten en vez de aplicarla sobre la definicioacuten del problema y su interpretacioacuten
El trabajo monoacutetono que se hacia anteriormente al uso de la computadora hace de importancia el dominio de los meacutetodos numeacutericos los cuales se deben llevar a cabo en combinacioacuten con las capacidades y potencialidades de la programacioacuten de computadoras para de esa forma resolver los problemas de ingenieriacutea mucho mas faacutecilmente y eficientemente
Objetivo de los meacutetodos de Runge-Kutta
El objetivo de los meacutetodos numeacutericos de runge-kutta es el anaacutelisis y solucioacuten de los problemas de valor inicial de ecuaciones diferenciales ordinarias (EDO) estos son una extensioacuten del meacutetodo de euler para resolver las (EDOrsquoS) pero con un orden de exactitud mas alto que este
Modelos matemaacuteticos
Los conocimientos cientiacuteficos se usan rutinariamente por los ingenieros en el disentildeo de elementos tales como maquinas circuitos eleacutectricos estructuras etc
Estos conocimientos son muy uacutetiles cuando se expresan en forma de un modelo matemaacutetico el cual se puede definir como una ecuacioacuten que expresa las caracteriacutesticas fundamentales de un sistema o proceso fiacutesico en teacuterminos matemaacuteticos siendo clasificados estos modelos desde simples relaciones algebraicas hasta grandes y complicados sistemas de ecuaciones diferenciales
Anaacutelisis de un modelo matemaacutetico
Un modelo matemaacutetico es una expresioacuten matemaacutetica como veremos en el siguiente ejemplo
Formula de la segunda ley de newton
F= ma donde F es la fuerza neta que actuacutea sobre el cuerpo
m es la masa del objeto
Utilizando esta ley vamos a determinar la velocidad de un paracaidista en caiacuteda libre
Para este caso puede crearse un nuevo modelo expresando la aceleracioacuten como la razoacuten de cambio de la velocidad con respecto al tiempo (dvdt)
Y sustituir en la ecuacioacuten de nueva forma
F=m(dvdt)
Asiacute la masa multiplicada por la razoacuten de cambio de la velocidad es igual a la suma de fuerzas que actuacutean 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 permaneceraacute a un nivel constante
Para un cuerpo que cae la fuerza la fuerza total esta compuesta por dos fuerzas contrarias la atraccioacuten 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 gravitacioacuten que equivale a 980 cm por segundo cuadrado
La resistencia del aire se puede formular como una aproximacioacuten 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 asiacute que combinando las ecuaciones anteriores
M(dvdt)= mg ndash cv o dividiendo cada lado entre m
esta es la ecuacioacuten de un cuerpo que cae a las fuerzas que actuacutean sobre el y es una ecuacioacuten diferencial porque esta escrita en teacuterminos de la razoacuten diferencial dvdt
Usando el calculo y resolviendo esta ecuacioacuten diferencial se puede llegar a la siguiente formula que expresa la velocidad del paracaidista en funcioacuten del tiempo
Ec(principal)
A continuacioacuten 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 demaacutes
Enunciado del problema de un paracaidista que cae
MEacuteTODO ANALIacuteTICO
Un paracaidista con una masa de 55500 g salta de un aeroplano apliquese la ecuacioacuten principal para calcular la velocidad antes de abrir el paracaiacutedas El cociente de arrastre c es aproximadamente igual a 10500 gs
Solucioacuten
Al sustituir los valores de los paraacutemetros en la ecuacioacuten 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 continuacioacuten
TABLA DE RESULTADOS
Tiempo en segundos
Velocidad en cms
0 0
2 16317
4 27495
6 35151
8 40396
10 43988
12 46449
Al infinito 51800
Nota La escala de la velocidad en la graacutefica es de 1=1000
Si lo has notado para sacar el resultado de tus caacutelculos de la tabla y graacutefica anteriores necesitas estar sustituyendo en la formula de v(t) esto hace el meacutetodo analiacutetico cansado y repetitivo pero es una solucioacuten analiacutetica exacta porque satisface la ecuacioacuten diferencial original
SOLUCIOacuteN NUMEacuteRICA
La aproximacioacuten de la razoacuten de cambio de la velocidad con respecto al tiempo puede representarse de la siguiente forma
Y llegamos a ordenar esta ecuacioacuten de la siguiente manera
la cual puede emplearse para extender el calculo asiacute que aplicandola directamente se diraacute
Nuevo valor = valor anterior + valor estimulado x incremento del tiempo
de v de v de la pendiente
Asiacute pues efectuando el mismo calculo para el problema se procede como sigue Sustituyendo en la ecuacioacuten nuevamente vemos que necesitamos aplicar las matemaacuteticas sin el auxilio de una computadora y nos vamos a tardar un buen rato
v(4)= 23308 cms
Y asiacute sucesivamente se continua con el calculo para obtener los valores en la siguiente tabla
Tiempo en seg
Vel En cms
0
2
4
6
8
10
12
al infinito
Puede verse por el resultado en la tabla anterior que la solucioacuten por un meacutetodo numeacuterico se aproxima bastante bien a la solucioacuten exacta pero debido al empleo de rectas para aproximar la funcioacuten que es continuamente curva existe discrepancia entre los resultados de la tabla del meacutetodo analiacutetico y la de este meacutetodo una manera de minimizar el error es utilizando intervalos menores de
tiempo en el muestreo de la ecuacioacuten y asiacute los segmentos de recta seguiraacuten mas de cerca a la solucioacuten verdadera
Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequentildeos haraacute poco practicas las soluciones numeacutericas es por eso que a continuacioacuten se presentara un meacutetodo con el cual se haraacute uso de un programa de computadora para la solucioacuten del problema
SOLUCIOacuteN POR EL MEacuteTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA
Los meacutetodos de euler son procedimientos para resolver EDOrsquoS de primer orden y se pueden programar en una computadora para no hacer los caacutelculos manualmente
Una desventaja importante de los meacutetodos 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 aplicacioacuten dentro de la programacioacuten por computadora para la solucioacuten del problemas
PROGRAMA PARA LA SOLUCIOacuteN DEL PROBLEMA DEL PARACAIDISTA
clear clf hold off t=0 n=0 v=0 C=105 M=555 g=98 h=1 t_rec(1)=t v_rec(1)=v while tlt=12 n=n+1 v=v+h(-CMvv+g) t=t+h v_rec(n+1)=v t_rec(n+1)=t end plot(t_recv_rec)
Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programacioacuten como fortran o C pero mas sencillo que estos puesto que cualquier variable puede contener nuacutemeros de cualquier tipo sin una declaracioacuten especial durante la programacioacuten
Por lo que resulta una poderosa herramienta para el anaacutelisis numeacuterico y matemaacutetico en la aplicacioacuten de los meacutetodos numeacutericos a la solucioacuten de los problemas de ingenieriacutea
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
512 Meacutetodo de Runge-Kutta
Introduccioacuten La computadora es la herramienta mas poderosa hasta ahora conocida para la solucioacuten de problemas en el campo de las ciencias exactas en este caso los meacutetodos numeacutericos como punto principal por sus aplicaciones en la ingenieriacutea
Los meacutetodos numeacutericos son teacutecnicas donde es posible resolver los problemas por medio de operaciones aritmeacuteticas estos meacutetodos implementan un buen numero de caacutelculos que son por demaacutes demasiado lentos si se hacen manualmente gastando mucha energiacutea en la teacutecnica misma de solucioacuten en vez de aplicarla sobre la definicioacuten del problema y su interpretacioacuten
El trabajo monoacutetono que se hacia anteriormente al uso de la computadora hace de importancia el dominio de los meacutetodos numeacutericos los cuales se deben llevar a cabo en combinacioacuten con las capacidades y potencialidades de la programacioacuten de computadoras para de esa forma resolver los problemas de ingenieriacutea mucho mas faacutecilmente y eficientemente
Objetivo de los meacutetodos de Runge-Kutta
El objetivo de los meacutetodos numeacutericos de runge-kutta es el anaacutelisis y solucioacuten de los problemas de valor inicial de ecuaciones diferenciales ordinarias (EDO) estos son una extensioacuten del meacutetodo de euler para resolver las (EDOrsquoS) pero con un orden de exactitud mas alto que este
Modelos matemaacuteticos
Los conocimientos cientiacuteficos se usan rutinariamente por los ingenieros en el disentildeo de elementos tales como maquinas circuitos eleacutectricos estructuras etc
Estos conocimientos son muy uacutetiles cuando se expresan en forma de un modelo matemaacutetico el cual se puede definir como una ecuacioacuten que expresa las caracteriacutesticas fundamentales de un sistema o proceso fiacutesico en teacuterminos matemaacuteticos siendo clasificados estos modelos desde simples relaciones algebraicas hasta grandes y complicados sistemas de ecuaciones diferenciales
Anaacutelisis de un modelo matemaacutetico
Un modelo matemaacutetico es una expresioacuten matemaacutetica como veremos en el siguiente ejemplo
Formula de la segunda ley de newton
F= ma donde F es la fuerza neta que actuacutea sobre el cuerpo
m es la masa del objeto
Utilizando esta ley vamos a determinar la velocidad de un paracaidista en caiacuteda libre
Para este caso puede crearse un nuevo modelo expresando la aceleracioacuten como la razoacuten de cambio de la velocidad con respecto al tiempo (dvdt)
Y sustituir en la ecuacioacuten de nueva forma
F=m(dvdt)
Asiacute la masa multiplicada por la razoacuten de cambio de la velocidad es igual a la suma de fuerzas que actuacutean 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 permaneceraacute a un nivel constante
Para un cuerpo que cae la fuerza la fuerza total esta compuesta por dos fuerzas contrarias la atraccioacuten 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 gravitacioacuten que equivale a 980 cm por segundo cuadrado
La resistencia del aire se puede formular como una aproximacioacuten 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 asiacute que combinando las ecuaciones anteriores
M(dvdt)= mg ndash cv o dividiendo cada lado entre m
esta es la ecuacioacuten de un cuerpo que cae a las fuerzas que actuacutean sobre el y es una ecuacioacuten diferencial porque esta escrita en teacuterminos de la razoacuten diferencial dvdt
Usando el calculo y resolviendo esta ecuacioacuten diferencial se puede llegar a la siguiente formula que expresa la velocidad del paracaidista en funcioacuten del tiempo
Ec(principal)
A continuacioacuten 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 demaacutes
Enunciado del problema de un paracaidista que cae
MEacuteTODO ANALIacuteTICO
Un paracaidista con una masa de 55500 g salta de un aeroplano apliquese la ecuacioacuten principal para calcular la velocidad antes de abrir el paracaiacutedas El cociente de arrastre c es aproximadamente igual a 10500 gs
Solucioacuten
Al sustituir los valores de los paraacutemetros en la ecuacioacuten 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 continuacioacuten
TABLA DE RESULTADOS
Tiempo en segundos
Velocidad en cms
0 0
2 16317
4 27495
6 35151
8 40396
10 43988
12 46449
Al infinito 51800
Nota La escala de la velocidad en la graacutefica es de 1=1000
Si lo has notado para sacar el resultado de tus caacutelculos de la tabla y graacutefica anteriores necesitas estar sustituyendo en la formula de v(t) esto hace el meacutetodo analiacutetico cansado y repetitivo pero es una solucioacuten analiacutetica exacta porque satisface la ecuacioacuten diferencial original
SOLUCIOacuteN NUMEacuteRICA
La aproximacioacuten de la razoacuten de cambio de la velocidad con respecto al tiempo puede representarse de la siguiente forma
Y llegamos a ordenar esta ecuacioacuten de la siguiente manera
la cual puede emplearse para extender el calculo asiacute que aplicandola directamente se diraacute
Nuevo valor = valor anterior + valor estimulado x incremento del tiempo
de v de v de la pendiente
Asiacute pues efectuando el mismo calculo para el problema se procede como sigue Sustituyendo en la ecuacioacuten nuevamente vemos que necesitamos aplicar las matemaacuteticas sin el auxilio de una computadora y nos vamos a tardar un buen rato
v(4)= 23308 cms
Y asiacute sucesivamente se continua con el calculo para obtener los valores en la siguiente tabla
Tiempo en seg
Vel En cms
0
2
4
6
8
10
12
al infinito
Puede verse por el resultado en la tabla anterior que la solucioacuten por un meacutetodo numeacuterico se aproxima bastante bien a la solucioacuten exacta pero debido al empleo de rectas para aproximar la funcioacuten que es continuamente curva existe discrepancia entre los resultados de la tabla del meacutetodo analiacutetico y la de este meacutetodo una manera de minimizar el error es utilizando intervalos menores de
tiempo en el muestreo de la ecuacioacuten y asiacute los segmentos de recta seguiraacuten mas de cerca a la solucioacuten verdadera
Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequentildeos haraacute poco practicas las soluciones numeacutericas es por eso que a continuacioacuten se presentara un meacutetodo con el cual se haraacute uso de un programa de computadora para la solucioacuten del problema
SOLUCIOacuteN POR EL MEacuteTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA
Los meacutetodos de euler son procedimientos para resolver EDOrsquoS de primer orden y se pueden programar en una computadora para no hacer los caacutelculos manualmente
Una desventaja importante de los meacutetodos 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 aplicacioacuten dentro de la programacioacuten por computadora para la solucioacuten del problemas
PROGRAMA PARA LA SOLUCIOacuteN DEL PROBLEMA DEL PARACAIDISTA
clear clf hold off t=0 n=0 v=0 C=105 M=555 g=98 h=1 t_rec(1)=t v_rec(1)=v while tlt=12 n=n+1 v=v+h(-CMvv+g) t=t+h v_rec(n+1)=v t_rec(n+1)=t end plot(t_recv_rec)
Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programacioacuten como fortran o C pero mas sencillo que estos puesto que cualquier variable puede contener nuacutemeros de cualquier tipo sin una declaracioacuten especial durante la programacioacuten
Por lo que resulta una poderosa herramienta para el anaacutelisis numeacuterico y matemaacutetico en la aplicacioacuten de los meacutetodos numeacutericos a la solucioacuten de los problemas de ingenieriacutea
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
m es la masa del objeto
Utilizando esta ley vamos a determinar la velocidad de un paracaidista en caiacuteda libre
Para este caso puede crearse un nuevo modelo expresando la aceleracioacuten como la razoacuten de cambio de la velocidad con respecto al tiempo (dvdt)
Y sustituir en la ecuacioacuten de nueva forma
F=m(dvdt)
Asiacute la masa multiplicada por la razoacuten de cambio de la velocidad es igual a la suma de fuerzas que actuacutean 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 permaneceraacute a un nivel constante
Para un cuerpo que cae la fuerza la fuerza total esta compuesta por dos fuerzas contrarias la atraccioacuten 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 gravitacioacuten que equivale a 980 cm por segundo cuadrado
La resistencia del aire se puede formular como una aproximacioacuten 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 asiacute que combinando las ecuaciones anteriores
M(dvdt)= mg ndash cv o dividiendo cada lado entre m
esta es la ecuacioacuten de un cuerpo que cae a las fuerzas que actuacutean sobre el y es una ecuacioacuten diferencial porque esta escrita en teacuterminos de la razoacuten diferencial dvdt
Usando el calculo y resolviendo esta ecuacioacuten diferencial se puede llegar a la siguiente formula que expresa la velocidad del paracaidista en funcioacuten del tiempo
Ec(principal)
A continuacioacuten 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 demaacutes
Enunciado del problema de un paracaidista que cae
MEacuteTODO ANALIacuteTICO
Un paracaidista con una masa de 55500 g salta de un aeroplano apliquese la ecuacioacuten principal para calcular la velocidad antes de abrir el paracaiacutedas El cociente de arrastre c es aproximadamente igual a 10500 gs
Solucioacuten
Al sustituir los valores de los paraacutemetros en la ecuacioacuten 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 continuacioacuten
TABLA DE RESULTADOS
Tiempo en segundos
Velocidad en cms
0 0
2 16317
4 27495
6 35151
8 40396
10 43988
12 46449
Al infinito 51800
Nota La escala de la velocidad en la graacutefica es de 1=1000
Si lo has notado para sacar el resultado de tus caacutelculos de la tabla y graacutefica anteriores necesitas estar sustituyendo en la formula de v(t) esto hace el meacutetodo analiacutetico cansado y repetitivo pero es una solucioacuten analiacutetica exacta porque satisface la ecuacioacuten diferencial original
SOLUCIOacuteN NUMEacuteRICA
La aproximacioacuten de la razoacuten de cambio de la velocidad con respecto al tiempo puede representarse de la siguiente forma
Y llegamos a ordenar esta ecuacioacuten de la siguiente manera
la cual puede emplearse para extender el calculo asiacute que aplicandola directamente se diraacute
Nuevo valor = valor anterior + valor estimulado x incremento del tiempo
de v de v de la pendiente
Asiacute pues efectuando el mismo calculo para el problema se procede como sigue Sustituyendo en la ecuacioacuten nuevamente vemos que necesitamos aplicar las matemaacuteticas sin el auxilio de una computadora y nos vamos a tardar un buen rato
v(4)= 23308 cms
Y asiacute sucesivamente se continua con el calculo para obtener los valores en la siguiente tabla
Tiempo en seg
Vel En cms
0
2
4
6
8
10
12
al infinito
Puede verse por el resultado en la tabla anterior que la solucioacuten por un meacutetodo numeacuterico se aproxima bastante bien a la solucioacuten exacta pero debido al empleo de rectas para aproximar la funcioacuten que es continuamente curva existe discrepancia entre los resultados de la tabla del meacutetodo analiacutetico y la de este meacutetodo una manera de minimizar el error es utilizando intervalos menores de
tiempo en el muestreo de la ecuacioacuten y asiacute los segmentos de recta seguiraacuten mas de cerca a la solucioacuten verdadera
Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequentildeos haraacute poco practicas las soluciones numeacutericas es por eso que a continuacioacuten se presentara un meacutetodo con el cual se haraacute uso de un programa de computadora para la solucioacuten del problema
SOLUCIOacuteN POR EL MEacuteTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA
Los meacutetodos de euler son procedimientos para resolver EDOrsquoS de primer orden y se pueden programar en una computadora para no hacer los caacutelculos manualmente
Una desventaja importante de los meacutetodos 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 aplicacioacuten dentro de la programacioacuten por computadora para la solucioacuten del problemas
PROGRAMA PARA LA SOLUCIOacuteN DEL PROBLEMA DEL PARACAIDISTA
clear clf hold off t=0 n=0 v=0 C=105 M=555 g=98 h=1 t_rec(1)=t v_rec(1)=v while tlt=12 n=n+1 v=v+h(-CMvv+g) t=t+h v_rec(n+1)=v t_rec(n+1)=t end plot(t_recv_rec)
Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programacioacuten como fortran o C pero mas sencillo que estos puesto que cualquier variable puede contener nuacutemeros de cualquier tipo sin una declaracioacuten especial durante la programacioacuten
Por lo que resulta una poderosa herramienta para el anaacutelisis numeacuterico y matemaacutetico en la aplicacioacuten de los meacutetodos numeacutericos a la solucioacuten de los problemas de ingenieriacutea
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
Usando el calculo y resolviendo esta ecuacioacuten diferencial se puede llegar a la siguiente formula que expresa la velocidad del paracaidista en funcioacuten del tiempo
Ec(principal)
A continuacioacuten 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 demaacutes
Enunciado del problema de un paracaidista que cae
MEacuteTODO ANALIacuteTICO
Un paracaidista con una masa de 55500 g salta de un aeroplano apliquese la ecuacioacuten principal para calcular la velocidad antes de abrir el paracaiacutedas El cociente de arrastre c es aproximadamente igual a 10500 gs
Solucioacuten
Al sustituir los valores de los paraacutemetros en la ecuacioacuten 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 continuacioacuten
TABLA DE RESULTADOS
Tiempo en segundos
Velocidad en cms
0 0
2 16317
4 27495
6 35151
8 40396
10 43988
12 46449
Al infinito 51800
Nota La escala de la velocidad en la graacutefica es de 1=1000
Si lo has notado para sacar el resultado de tus caacutelculos de la tabla y graacutefica anteriores necesitas estar sustituyendo en la formula de v(t) esto hace el meacutetodo analiacutetico cansado y repetitivo pero es una solucioacuten analiacutetica exacta porque satisface la ecuacioacuten diferencial original
SOLUCIOacuteN NUMEacuteRICA
La aproximacioacuten de la razoacuten de cambio de la velocidad con respecto al tiempo puede representarse de la siguiente forma
Y llegamos a ordenar esta ecuacioacuten de la siguiente manera
la cual puede emplearse para extender el calculo asiacute que aplicandola directamente se diraacute
Nuevo valor = valor anterior + valor estimulado x incremento del tiempo
de v de v de la pendiente
Asiacute pues efectuando el mismo calculo para el problema se procede como sigue Sustituyendo en la ecuacioacuten nuevamente vemos que necesitamos aplicar las matemaacuteticas sin el auxilio de una computadora y nos vamos a tardar un buen rato
v(4)= 23308 cms
Y asiacute sucesivamente se continua con el calculo para obtener los valores en la siguiente tabla
Tiempo en seg
Vel En cms
0
2
4
6
8
10
12
al infinito
Puede verse por el resultado en la tabla anterior que la solucioacuten por un meacutetodo numeacuterico se aproxima bastante bien a la solucioacuten exacta pero debido al empleo de rectas para aproximar la funcioacuten que es continuamente curva existe discrepancia entre los resultados de la tabla del meacutetodo analiacutetico y la de este meacutetodo una manera de minimizar el error es utilizando intervalos menores de
tiempo en el muestreo de la ecuacioacuten y asiacute los segmentos de recta seguiraacuten mas de cerca a la solucioacuten verdadera
Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequentildeos haraacute poco practicas las soluciones numeacutericas es por eso que a continuacioacuten se presentara un meacutetodo con el cual se haraacute uso de un programa de computadora para la solucioacuten del problema
SOLUCIOacuteN POR EL MEacuteTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA
Los meacutetodos de euler son procedimientos para resolver EDOrsquoS de primer orden y se pueden programar en una computadora para no hacer los caacutelculos manualmente
Una desventaja importante de los meacutetodos 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 aplicacioacuten dentro de la programacioacuten por computadora para la solucioacuten del problemas
PROGRAMA PARA LA SOLUCIOacuteN DEL PROBLEMA DEL PARACAIDISTA
clear clf hold off t=0 n=0 v=0 C=105 M=555 g=98 h=1 t_rec(1)=t v_rec(1)=v while tlt=12 n=n+1 v=v+h(-CMvv+g) t=t+h v_rec(n+1)=v t_rec(n+1)=t end plot(t_recv_rec)
Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programacioacuten como fortran o C pero mas sencillo que estos puesto que cualquier variable puede contener nuacutemeros de cualquier tipo sin una declaracioacuten especial durante la programacioacuten
Por lo que resulta una poderosa herramienta para el anaacutelisis numeacuterico y matemaacutetico en la aplicacioacuten de los meacutetodos numeacutericos a la solucioacuten de los problemas de ingenieriacutea
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
8 40396
10 43988
12 46449
Al infinito 51800
Nota La escala de la velocidad en la graacutefica es de 1=1000
Si lo has notado para sacar el resultado de tus caacutelculos de la tabla y graacutefica anteriores necesitas estar sustituyendo en la formula de v(t) esto hace el meacutetodo analiacutetico cansado y repetitivo pero es una solucioacuten analiacutetica exacta porque satisface la ecuacioacuten diferencial original
SOLUCIOacuteN NUMEacuteRICA
La aproximacioacuten de la razoacuten de cambio de la velocidad con respecto al tiempo puede representarse de la siguiente forma
Y llegamos a ordenar esta ecuacioacuten de la siguiente manera
la cual puede emplearse para extender el calculo asiacute que aplicandola directamente se diraacute
Nuevo valor = valor anterior + valor estimulado x incremento del tiempo
de v de v de la pendiente
Asiacute pues efectuando el mismo calculo para el problema se procede como sigue Sustituyendo en la ecuacioacuten nuevamente vemos que necesitamos aplicar las matemaacuteticas sin el auxilio de una computadora y nos vamos a tardar un buen rato
v(4)= 23308 cms
Y asiacute sucesivamente se continua con el calculo para obtener los valores en la siguiente tabla
Tiempo en seg
Vel En cms
0
2
4
6
8
10
12
al infinito
Puede verse por el resultado en la tabla anterior que la solucioacuten por un meacutetodo numeacuterico se aproxima bastante bien a la solucioacuten exacta pero debido al empleo de rectas para aproximar la funcioacuten que es continuamente curva existe discrepancia entre los resultados de la tabla del meacutetodo analiacutetico y la de este meacutetodo una manera de minimizar el error es utilizando intervalos menores de
tiempo en el muestreo de la ecuacioacuten y asiacute los segmentos de recta seguiraacuten mas de cerca a la solucioacuten verdadera
Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequentildeos haraacute poco practicas las soluciones numeacutericas es por eso que a continuacioacuten se presentara un meacutetodo con el cual se haraacute uso de un programa de computadora para la solucioacuten del problema
SOLUCIOacuteN POR EL MEacuteTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA
Los meacutetodos de euler son procedimientos para resolver EDOrsquoS de primer orden y se pueden programar en una computadora para no hacer los caacutelculos manualmente
Una desventaja importante de los meacutetodos 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 aplicacioacuten dentro de la programacioacuten por computadora para la solucioacuten del problemas
PROGRAMA PARA LA SOLUCIOacuteN DEL PROBLEMA DEL PARACAIDISTA
clear clf hold off t=0 n=0 v=0 C=105 M=555 g=98 h=1 t_rec(1)=t v_rec(1)=v while tlt=12 n=n+1 v=v+h(-CMvv+g) t=t+h v_rec(n+1)=v t_rec(n+1)=t end plot(t_recv_rec)
Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programacioacuten como fortran o C pero mas sencillo que estos puesto que cualquier variable puede contener nuacutemeros de cualquier tipo sin una declaracioacuten especial durante la programacioacuten
Por lo que resulta una poderosa herramienta para el anaacutelisis numeacuterico y matemaacutetico en la aplicacioacuten de los meacutetodos numeacutericos a la solucioacuten de los problemas de ingenieriacutea
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
de v de v de la pendiente
Asiacute pues efectuando el mismo calculo para el problema se procede como sigue Sustituyendo en la ecuacioacuten nuevamente vemos que necesitamos aplicar las matemaacuteticas sin el auxilio de una computadora y nos vamos a tardar un buen rato
v(4)= 23308 cms
Y asiacute sucesivamente se continua con el calculo para obtener los valores en la siguiente tabla
Tiempo en seg
Vel En cms
0
2
4
6
8
10
12
al infinito
Puede verse por el resultado en la tabla anterior que la solucioacuten por un meacutetodo numeacuterico se aproxima bastante bien a la solucioacuten exacta pero debido al empleo de rectas para aproximar la funcioacuten que es continuamente curva existe discrepancia entre los resultados de la tabla del meacutetodo analiacutetico y la de este meacutetodo una manera de minimizar el error es utilizando intervalos menores de
tiempo en el muestreo de la ecuacioacuten y asiacute los segmentos de recta seguiraacuten mas de cerca a la solucioacuten verdadera
Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequentildeos haraacute poco practicas las soluciones numeacutericas es por eso que a continuacioacuten se presentara un meacutetodo con el cual se haraacute uso de un programa de computadora para la solucioacuten del problema
SOLUCIOacuteN POR EL MEacuteTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA
Los meacutetodos de euler son procedimientos para resolver EDOrsquoS de primer orden y se pueden programar en una computadora para no hacer los caacutelculos manualmente
Una desventaja importante de los meacutetodos 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 aplicacioacuten dentro de la programacioacuten por computadora para la solucioacuten del problemas
PROGRAMA PARA LA SOLUCIOacuteN DEL PROBLEMA DEL PARACAIDISTA
clear clf hold off t=0 n=0 v=0 C=105 M=555 g=98 h=1 t_rec(1)=t v_rec(1)=v while tlt=12 n=n+1 v=v+h(-CMvv+g) t=t+h v_rec(n+1)=v t_rec(n+1)=t end plot(t_recv_rec)
Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programacioacuten como fortran o C pero mas sencillo que estos puesto que cualquier variable puede contener nuacutemeros de cualquier tipo sin una declaracioacuten especial durante la programacioacuten
Por lo que resulta una poderosa herramienta para el anaacutelisis numeacuterico y matemaacutetico en la aplicacioacuten de los meacutetodos numeacutericos a la solucioacuten de los problemas de ingenieriacutea
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
tiempo en el muestreo de la ecuacioacuten y asiacute los segmentos de recta seguiraacuten mas de cerca a la solucioacuten verdadera
Calculando manualmente el esfuerzo al usar incrementos cada vez mas pequentildeos haraacute poco practicas las soluciones numeacutericas es por eso que a continuacioacuten se presentara un meacutetodo con el cual se haraacute uso de un programa de computadora para la solucioacuten del problema
SOLUCIOacuteN POR EL MEacuteTODO DE EULER DEL PROBLEMA DEL PARACAIDISTA
Los meacutetodos de euler son procedimientos para resolver EDOrsquoS de primer orden y se pueden programar en una computadora para no hacer los caacutelculos manualmente
Una desventaja importante de los meacutetodos 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 aplicacioacuten dentro de la programacioacuten por computadora para la solucioacuten del problemas
PROGRAMA PARA LA SOLUCIOacuteN DEL PROBLEMA DEL PARACAIDISTA
clear clf hold off t=0 n=0 v=0 C=105 M=555 g=98 h=1 t_rec(1)=t v_rec(1)=v while tlt=12 n=n+1 v=v+h(-CMvv+g) t=t+h v_rec(n+1)=v t_rec(n+1)=t end plot(t_recv_rec)
Este programa anterior fue realizado en matlab (Laboratorio de matrices) el cual puede considerarse como un lenguaje de programacioacuten como fortran o C pero mas sencillo que estos puesto que cualquier variable puede contener nuacutemeros de cualquier tipo sin una declaracioacuten especial durante la programacioacuten
Por lo que resulta una poderosa herramienta para el anaacutelisis numeacuterico y matemaacutetico en la aplicacioacuten de los meacutetodos numeacutericos a la solucioacuten de los problemas de ingenieriacutea
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
La siguiente graacutefica muestra el resultado del programa anterior (tambieacuten realizada por matlab)
Velocidad cms Tiempo en seg
Podemos observar con los resultados anteriores que la solucioacuten por medio de la programacioacuten por computadora se puede hacer que se aproxime a la solucioacuten exacta de la ecuacioacuten diferencial tanto como se quiera dependiendo del meacutetodo que se utilice
En los meacutetodos de RUNGE-KUTTA el meacutetodo de exactitud se incrementa mediante el empleo de un meacutetodo de integracioacuten numeacuterica 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 ecuacioacuten diferencial ordinaria siguiente
Si queremos calcular
Aplicamos una regla trapezoidal al miembro derecho la cual resulta de dividir el aacuterea de una curva en n intervalos iguales de longitud h evaluando el aacuterea de cada una de las fajas aproximando el arco y=f(x) con la cuerda correspondiente a cada fajita si h es lo suficientemente pequentildea entonces se puede aproximar como el aacuterea de un trapecio bajo la curva es decir
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
El aacuterea total es la suma de las aacutereas parciales
Entonces la forma trapezoidal en donde el aacuterea 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 meacutetodos de RUNGE-KUTTA para poder aplicar las ecuaciones a problemas concretos
A continuacioacuten veremos algunos ejemplos de los meacutetodos de runge a problemas de circuitos eleacutectricos y su solucioacuten mas raacutepida 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 eleacutectrico cambia significativamente dependiendo de los valores de los componentes empleados asiacute en el circuito que se muestra a continuacioacuten 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 ecuacioacuten diferencial
I(0) = 0
Se necesita encontrar el valor de la corriente para 0lt t lt 002 s mediante un meacutetodo de RUNGE-KUTTA de segundo orden con h=0001
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
donde
Mostramos entonces los caacutelculos para los dos primeros pasos
t = 0
t= 00001
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 ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
tardado resolverlo manualmente de ahiacute que lo que falta para calcular se haraacute en el programa de matlab
NOTA Los programas pueden estar sujetos a errores al escribirse asiacute que es recomendable estudiar matlab con mayor detenimiento para poder detectarlos mas faacutecilmente
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=RL EL=EL while t(n)lt02 K1=hfn10_9(y(n)RLEL) K2=hfn10_9(y(n)+K1RLEL) y(n+1)=y(n)+5(K1+K2) t(n+1)=n+h n=n+1 end plot(ty) xlabel(t = seg ) ylabel(I(A)) fn10_9 function f=fn10_9(IRLEL) f=-RLI+EL
Ahora vamos a un segundo ejemplo de aplicacioacuten a circuitos y veamos la siguiente figura
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
Los valores de los elementos que se encuentran en la figura anterior son hipoteacuteticos es solo para representar el circuito el cual puede definirse por medio de las siguientes ecuaciones diferenciales con su consiguiente solucioacuten por medio de un listado de programa Para la primera inductancia resistencia y capacitor
con
donde e(t) = 0 excepto e(t)=1 cuando 0lttlt001s 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 solucioacuten se determina en el siguiente programa realizado en matlab
clearclg subplot(221) for k=14 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=20C=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=[-RaLa RaLa -1(LaC) RaLb-(Ra+Rb)Lb 1(LbC)
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
1C-1C0] S=[000] X=[000] h=00005 for n=1101 t=(n-1)h S=[sin(t600)exp(-t600)La00] S=[cos(t600)La00] S=[1La00] if tgt001 S=[000]end k1=h(MX+S) k2=h(M(X+k12)+S) k3=h(M(X+k22)+S) k4=h(M(X+k3)+S) X=X+(K1+K22+k32+k4)6 X_r(n)=X t_r(n)=t end
plot(t_rX_r(12)t_rX_r(3)) xlabel(t)ylabel(i1i2q) L=length(t_r) text(t_r(L10)X_r(1L10)i1) text(t_r(L2)X_r(2L2)i2) text(t_r(L8)X_r(3L8)q) if k==1title((A))end if k==2title((B))end if k==3title((C))end if k==4title((D))end end
Hay que tener en cuenta que de alguna u otra manera los meacutetodos numeacutericos estan sujetos a ciertos errores al hacer los caacutelculos 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 pequentildeos para la aplicacioacuten de las ecuaciones y obtener los resultados mas raacutepido y de manera mas exacta
No hay limite para las aplicaciones que se pueden realizar en programacioacuten de computadoras para resolver los problemas de meacutetodos numeacutericos asiacute que ya para terminar a continuacioacuten se daraacuten algunos otros ejemplos en matlab para algunos meacutetodos mencionados como claacutesicos en algunos libros
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
52 Meacutetodo de pasos muacuteltiples
Meacutetodos de Paso MultipleEL valor de xi+1se realiza en funcioacuten de los valores de xixi-1 hellipxi-p
El paso de integracioacuten h puede ser fijo o variable En caso de variabilidad se utilizaraacute el error de integracioacuten estimado como diferencia en la evaluacioacuten de xi+1con dos meacutetodos diferentesMetodos de paso multipleSe considera el problema de valores iniciales (PVI)8lty0(x) = f(x y(x)) x 2 [a b]y(a) = y0 dadoel que supondremos tiene solucion unica y [a b] 1048576 R la cual es acotada y depende continuamente de los datos f e y0Dada una partici acuteon del intervalo [a b]a = x0 lt x1 lt _ _ _ lt xN = blos metodos que hemos visto hasta aqui solo usan la informacion del valor y de la solucion
calculada en xi para obtener yi+1 Por eso se denominan metodos de paso simple
Parece razonable pensar que tambien podrian utilizarse los valores yi1048576k yiobtenidos en los nodos xi1048576k xi para k gt 0Para ello si integramos y0(x) = f(x y(x)) en el intervalo [xi xi+1] se tieneZ xi+1xiy0(x) dx = Z xi+1xif(x y(x)) dxy por lo tantoy(xi+1) = y(xi) + Z xi+1xif(x y(x)) dxLa integral puede calcularse aproximando el integrando mediante el polinomio de
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias
interpolacion de f(x y(x)) en los puntos xi1048576k xi Sin embargo como los valoresexactos de y(xi1048576k) y(xi) no se conocen no podemos utilizarlos en la evaluacionde f(x y(x)) En su lugar podemos utilizar los valores de la soluci acuteon calculadayi1048576k yi
Los metodos que asi se obtienen se denominan metodos de paso multiple
53 Sistemas de ecuaciones diferenciales ordinarias