Cálculo Numérico y Computación
Universidad Nacional de Cuyo Facultad de Ingeniería
CÁLCULO NUMÉRICO Y COMPUTACIÓN.
MATERIAL PARA CLASES DE TEORÍA. El Material para Clases de Teoría, es el conjunto de diapositivas con el que se desarrollarán las clases de Teoría. Está ordenado según el cronograma a seguir en e presente año lectivo. En general No incluye las demostraciones a desarrollar. Pero si tiene muchos de los gráficos y tablas con las que se desarrollan las clases de teoría. En general las tablas están asociadas a ejercicios o conceptos de síntesis; y suelen estar en blanco para sean completadas durante las clases.
Los temas desarrollados son los siguientes: Algoritmia. Conceptos Básicos Solución Numérica de Raíces de Ecuaciones No Lineales Sistemas de Ecuaciones Lineales Valores y Vectores Propios Interpolación y Aproximación de funciones discretas Integración Numérica Derivación Numérica, con aplicación a la solución de Ecuaciones Diferenciales con Valores de Contorno Solución Numérica de Ecuaciones Diferenciales Ordinarias con Valores Iniciales El ordenamiento elegido está según el cronograma del año 2016.
Dr. Ing. Anibal Mirasso
Página 1 de 2 Año 2016
Cálculo Numérico y Computación
Dr. Ing. Anibal Mirasso
Universidad Nacional de Cuyo Facultad de Ingeniería
Página 2 de 2 Año 2016
Facultad de Ingeniería, Universidad Nacional de Cuyo
• Definición de Algoritmo • Ejemplo • Pseudo código, Sintaxis de MATLAB • Variables, definición y tipos • Estructuras Algorítmicas Secuencial Decisión Variar Mientras • Ejercicios • Tips • Subprogramas en MATLAB
Facultad de Ingeniería, Universidad Nacional de Cuyo
Es la descripción de un procedimiento en forma ordenada y sistemática en un número finito de pasos, líneas o comandos. Se puede hacer mediante Pseudo código; Diagrama de Flujo; Diagrama de Bloques o Código de Programación Ejemplo: Algoritmo para calcular el promedio de dos números
Pseudo código
Elementos
Programa prom_num
Nombre
Reales (xa; xb; prom)
Declaración de variables
Escribir “Ingrese el primer número” Leer xa
Ingreso de datos
Escribir “Ingrese el segundo número” Leer xb prom
(xa+xb)/2
Escribir “el Promedio es: “ Escribir prom Fin
Sintaxis de MATLAB function prom_num % Reales (xa; xb; prom) disp('Ingrese el número 1:'); xa = input(''); disp('Ingrese el número 2:'); xb = input(''); prom = (xa+xb)/2; disp('el Promedio es:'); disp(prom);
Proceso Entrega de resultados Final
end
Facultad de Ingeniería, Universidad Nacional de Cuyo
Variable es una dirección, alojada en la memoria de la computadora, para la identificación del CPU es una especie de recipiente en cuyo interior podemos colocar cierta información se identifica con un nombre representativo
altura a
Variables simples solo puede guardar un solo valor y son las variables que ya hemos vistos: altura; base, xa, xb, etc Variables Dimensionadas son las variables que guardan información referida a un mismo dato y que por la cantidad de datos es necesario dimensionarlas. Tipicamente Vectores y Matrices. Vectores. Se dimensionan con el máximo número de componentes en la forma A(N).Se identifica uno con A(j). Nombre de la variable dimensionada
Casillero donde van los distintos valores de la variable “A”
A(j) A
A(N)
Matrices. Se dimensionan con el máximo número de componentes en la forma A(Nf,Nc).Se identifica uno con A(i,j). k=1: n3 j=1: Nc
j=1: Nc B(2,5) del B(Nf,Nc)
i varía de 1 a Nf
i varía de 1 a Nf
C(2,5,1) del C(Nf,Nc;n3)
Facultad de Ingeniería, Universidad Nacional de Cuyo
Variable pueden ser clasificadas según el TIPO de información que almacenan: Reales
Enteras
Lógicas
Caracter
Jerarquía de las operaciones con las distintas variables Prioridad
Operador
Resultado
1 2 3
^ */ +-+
Potencia Producto - Cociente Suma – Resta – Concatenación de caracteres
4 5 6
= ≠ < ≤ > ≥ .NOT. .AND.
7
.OR.
Relación Negación Conjunción [Y] lógico Disyunción [O] lógico
Numérico Numérico Suma y resta el resultado es numérico. Concatenación el resultado es Carácter Lógico Lógico Lógico Lógico
Facultad de Ingeniería, Universidad Nacional de Cuyo
! !
!
Es un conjunto finito de líneas, con principio y fin Ejemplo: Algoritmo para calcular las raíces de la ecuación de segundo grado
Pseudo código
Elementos
Sintaxis de MATLAB
Programa raices
Nombre
function raices
Reales (a; b; c; r1; r2)
Declaración de variables
% Reales (a; b; c; r1; r2) disp('Ingrese Coef Cuadrático:'); a = input('');
Escribir “Ingrese Coef Cuadrático” Leer a
Ingreso de datos
Escribir “Ingrese Coef Lineal”
b = input('Ingrese Coef Lineal:'); c = input('Ingrese Coef Indep:');
Leer b Escribir “Ingrese Coef Indep”, Leer c Discrim
(b^2-4*a*c)
r1
(-b+Discrim^0.5)/(2*a)
r2
(-b-Discrim^0.5)/(2*a)
Escribir “Raices son: “, r1, r2 Fin
Discriminante puede ser Negativo!!
Proceso
Discrim = (b^2-4*a*c); r1 = (-b+Discrim^0.5)/(2*a); r2 = (-b-Discrim^0.5)/(2*a);
Entrega de resultados Final
disp('Raices son:');disp(r1);disp(r2);
end
Facultad de Ingeniería, Universidad Nacional de Cuyo
! ! Es un conjunto finito de líneas, con una bifurcación según la respuesta lógica a una pregunta Ejemplo: Algoritmo para calcular las raíces de la ecuación de segundo grado
Pseudo código
Elementos
Programa raices
Nombre
Reales (a; b; c; r1; r2)
Declaración de variables
Escribir “Ingrese Coef Cuadrático”, Leer a Escribir “Ingrese Coef Lineal”, Leer b Escribir “Ingrese Coef Indep”, Leer c
Discrim
(b^2-4*a*c)
Ingreso de datos
Sintaxis de MATLAB function raices % Reales (a; b; c; r1; r2) a = input('Ingrese Coef Cuadrático:'); b = input('Ingrese Coef Lineal:'); c = input('Ingrese Coef Indep');
Discrim = (b^2-4*a*c); If (Discrim > 0) r1 = (-b+Discrim^0.5)/(2*a); r2 = (-b-Discrim^0.5)/(2*a);
Proceso
SI (Discrim > 0) entonces r1 (-b+Discrim^0.5)/(2*a) r2 (-b-Discrim^0.5)/(2*a) Escribir “Raices son: “, r1, r2
y Entrega de resultados
disp('Raices son:');disp(r1);disp(r2);
else
Sino
disp('Raices imaginarias:');
Escribir “Raíces imaginarias“
end
Fin del SI Fin
Final
end
Facultad de Ingeniería, Universidad Nacional de Cuyo
! ! Es un conjunto finito de líneas, que se repite un número definido y conocido de veces Ejemplo: Algoritmo para leer y escribir las componentes de un vector de dimensión 10
Pseudo código
Elementos
Programa leer_vec
Nombre
Reales (vec(10))
Declaración de variables
Enteros i, j DOFOR i =1 HASTA 10 PASO 1 Escribir “Ingrese componente”, Leer vec(i)
Ingreso de datos
ENDDO DOFOR j =1 HASTA 10 PASO 1 Escribir “La componente es”, Escribir vec(j)
function leer_vec % Reales (vec(10)) % Enteros i, j disp(lectura de las componentes'); for i=1:10 disp ('Ingrese Componente:'); vec(i) = input('');
end Proceso y Entrega de resultados
ENDDO Fin
Sintaxis de MATLAB
Final
disp('Escritura de las componentes'); for j=1:10 disp ('La componente es:'); disp (vec(j));
end end
Facultad de Ingeniería, Universidad Nacional de Cuyo
! ! Es un conjunto finito de líneas, que se repite un número indefinido y desconocido de veces Ejemplo: Algoritmo para que:”mientras la función f(x)=(3x2-12) sea distinto de cero, lea un valor de abscisa x0, calcule la función f(x0) en ese valor x0, y sólo entregue un valor de x0 cuando la función f(x) sea cero”.
Pseudo código
Elementos
Programa leer_vec
Nombre
Reales (x0, f)
Declaración de variables
Enteros k k 1 Escribir “Ingrese x0”, Leer x0 f 3*x0^2-12 DOWHILE k k+1 Escribir “Ingrese x0”, Leer x0 f 3*x0^2-12 ENDDO
¿Cuál es la utilidad de k?
function leer_vec %Real (f, x0) %Integer k
k=1; x0=input ('íngrese abscisa'); f= 3*x0^2-12;
Ingreso de datos Proceso y Entrega de resultados
while
(abs(f) > 0)
k =k + 1; x0=input ('íngrese abscisa'); f= 3*x0^2-12 ;
end disp('la abscisa es: '),x0
Escribir “la abscisa que anula f es: “, x0 Fin
Sintaxis de MATLAB
Final
end
Modificar para tener la “historia” de los x0
Facultad de Ingeniería, Universidad Nacional de Cuyo
! ! Es un conjunto finito de líneas, que se repite un número indefinido y desconocido de veces Ejemplo: Algoritmo para que:”mientras la función f(x)=(3x2-12) sea distinto de cero, lea un valor de abscisa x0, calcule la función f(x0) en ese valor x0.Y sólo entregue la historia de x0 cuando la función f(x) es cero”. Pseudo código Elementos Sintaxis de MATLAB Programa cer_de_f
Reales (x0, f) Enteros k k 1 Escribir “Ingrese x0”, Leer x0 f 3*x0^2-12
Nombre
Declaración de variables
Ingreso de datos
A(k,1)=x0; A(k,2)=f;
DOWHILE k k+1 Escribir “Ingrese x0”, Leer x0 f 3*x0^2-12 A(k,1)=x0;
Proceso y Entrega de resultados
function cer_de_f
%Real (A(2,1000), x0, f) %Integer k k=1; x0=input ('íngrese abscisa'); f= 3*x0^2-12; A(k,1)=x0; A(k,2)=f;
while
(abs(f) > 0) k =k + 1; x0=input ('íngrese abscisa'); f= 3*x0^2-12 ; A(k,1)=x0; A(k,2)=f;
A(k,2)=f;
end
ENDDO Escribir “la historio de : “, x0
disp('la historia de x0 y f es: ') for i=1:k disp(' '), A(k,1), A(k,2)
for i=1:k Escribir, A(k,1), A(k,2) end
Fin
end Final
end
Modificar para controlar que no se supere la dimensión 1000 de la matriz A y definir “banda de cero”
Facultad de Ingeniería, Universidad Nacional de Cuyo
! ! Ejercicio Algoritmo para que:”mientras la función f(x)=(3x2-12) sea distinto de cero, lea un valor de abscisa x0, calcule la función f(x0) en ese valor x0.Y sólo entregue la historia de x0 cuando la función f(x) es cero”. No se debe superar las 1234 iteraciones. function cer_de_f
%Real (A(2,1000), x0, f, tol) %Integer k, i tol=1e-12; max= ; k=1; x=input ('íngrese abscisa'); f= 3*x^2-12; A(1,k)=x; A(2,k)=f;
while (comparación o variable lógica) Línea 1 Línea 2 ……… ………
end
Bloque de líneas que se “ejecuta” sólo si la comparación o variable lógica es VERDADERA. Debe iniciarse antes del while Debe cambiar dentro del Bloque
% .
(abs(f) > tol ……. k < max) k =k + 1; x=input ('íngrese componente'); f= 3*x^2-12 ; A(1,k)=x; A(2,k)=f; .
disp('la historia de x0 y f es: ') for i=1:k disp(' '), A(i,1), A(i,2)
end end
Operadores Lógicos Mayor, Menor,……..….en MATLAB……>…< Igual, Distinto……..….en MATLAB…....==……~= OR…………………….en MATLAB……| AND………..................en MATLAB…..& dan como resultado Verdadero o Falso
Facultad de Ingeniería, Universidad Nacional de Cuyo
!
"
Los subprogramas son un subproceso constituido por un bloque de órdenes o comandos que realizan dicho subproceso y se los identifica con un Nombre. Se los llama por su Nombre desde un proceso principal Pueden necesitar datos de ingreso, y en general entregan un resultado Los datos de ingreso y resultados pueden intercambiarse con el proceso principal mediante argumentos En MATLAB se los identifica como “function”. Programa Principal Es un conjunto de líneas de órdenes, y en alguna de ellas se invoca una function mediante el Nombre de la misma function --------------------------------------------------% ingreso de datos [x1, y1, x2, y2]= leer_puntos; % sin pasaje de argumentos de entrada % con multiples argumentos de salida --------------------------------------------------------% Cálculo del Punto Medio ym = Punto_medio(x1, y1, x2, y2); % con pasaje de argumentos de entrada % con un único argumento de salida -----------------------------
NOTAR Variables Globales y Locales
Subprograma las function puede estar en el archivo principal o en otro archivo que debe llamarse con el nombre de la function function [a, b, c, a=input('Ingrese b=input('Ingrese c=input('Ingrese d=input('Ingrese end
d]= leer_puntos x1'); y1'); x2'); y2')
function zm = Punto_medio (za, zb, zc , zd) xx=(za+zc)/2; zm=(zb+zd)/2; end
Facultad de Ingeniería, Universidad Nacional de Cuyo
# Ejercicio 1 Desarrollar un algoritmo que lea un vector de 100 componentes e imprima su módulo (norma cuadrática) Ejercicio 2 Desarrollar un algoritmo que lea un vector de 45 componentes e imprima su norma infinita Ejercicio 3 Desarrollar un algoritmo que lea un vector de 45 componentes e imprima su versor, sólo si su norma infnito es no nula. De lo contario exprese que el vector es nulo
• No hay obviedades • No debe haber ambigüedades • Cada orden deber ser precisa • Ordene cosas simples, que sea sólo una tarea • Deje la optimización de cantidad de órdenes para una segunda etapa • Proponga el algoritmo y “ejecútelo”. Sea estricto en ello. • Si aparecen errores, corrija y vuelva a ejecutar. Itere hasta que esté conforme.
Facultad de Ingeniería, Universidad Nacional de Cuyo
RAÍCES DE ECUACIONES NO LINEALES
Planteo de Problema Estrategias Iterativas en general Análisis de función Acercamiento o Inicialización Recurrencia Control de Detención Actualización Métodos de Intervalo Bisección Regula Falsi Métodos Abiertos Método de la Secante Método de Newton Raphson Métodos de Punto Fijo Ejemplos Resumen
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
1 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
RAÍCES DE ECUACIONES NO LINEALES Se busca una aproximación de la abscisa que hace nulo el valor de una función, con una precisión deseada. Dicha abscisa se denomina raíz de la ecuación no lineal. Se obtiene de imponer igual a cero la ordenada (imagen) de la función. Supuestos
La función y = f(x): R ĺ R, es no singular, continua, conocida en forma analítica y tiene al menos una raíz
Y(x) y=f(x)
X X= r Ejemplo 1 Dada la función La ecuación no lineal Tiene como raíces r1, 2
f(x)=0 La ecuación no lineal es La abscisa X=r es Raíz si verifica la ecuación no lineal. f(r)=0 Es decir si:
y ax bxc 0 a x2 b x c 2
( b r b 2 4 a c ) /( 2 a )
Ejemplo 2 Dada la función La ecuación no lineal
y x tan( x ) c 0 x tan( x ) c
Tiene como raíces!!!!!!!.???????
Objetivo
Calcular una aproximación de la raíz mediante Métodos Iterativos. Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
2 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
PROCEDIMIENTO GENERAL Paso Inicial Realizar un análisis de la función: determinar singularidades, asíntotas, etc. Se busca toda la información posible a los efectos de elegir adecuadamente las variables iniciales de los métodos iterativos Paso Acercamiento o Inicialización Se debe elegir adecuadamente las variables iniciales de los métodos iterativos buscando toda la información posible a tal efecto. Se trata de encontrar un intervalo en el eje X donde exista al menos una raíz de la ecuación no lineal. Se busca un valor de abscisa cercano a la raíz. Paso Aproximación mediante Métodos Iterativos Métodos de Intervalos o Cerrados Métodos Abiertos Método de Bisección
Método de la Secante o de Newton Lagrange
Método de Regula Falsi
Método de Newton Métodos de Puntos Fijos
Siempre requieren de alguna condición de inicialización
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
3 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
Paso Acercamiento o Inicialización Se trata de encontrar un intervalo en el eje X donde exista al menos una raíz de la ecuación no lineal. Y(x)
Y(x) Y=f(x)
Y=f(x)
f(bk)
f(bk) f(ak)
X
ak f(ak)
bk
X= r
ak
Y=f(x) f(bk)
ak f(ak)
Dr. Ing. A.Mirasso
X= r1
X= r3
X= r2 X bk
f (a k ) f (bk ) 0
Y(x)
X= r2
X= r1
X
existe al menos una raìz
bk
Raíces de Ecuaciones No Lineales
Año 2014
4 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
Paso Aproximación mediante Métodos Iterativos INICIALIZACIÓN Se deben definir los contenidos de las variables de modo que se cumplan las condiciones de Inicialización del método HACER MIENTRAS “No Hay Solución” es Verdadero DOWHILE (NHS) RECUERRENCIA Se debe evaluar la nueva solución aproximada correspondiente a la nueva iteración o ciclo CONTROL DE DETENCIÓN Si alguna MEDIDA del ERROR es adecuada entonces se ha logrado la solución buscada. Se debe Cambiar NHS a falso SI ( Valor Absoluto de NHS es FALSO FINSI
f(rk+1) < Tolerancia)
ACTUALIZACIÓN DE VARIABLES Se reasignan las variable de modo que se cumplan con las condiciones de Inicalización FIN DEL MIENTRA
MÉTODOS DE INTERVALOS O CERRADOS Método de Bisección Método de Regula Falsi Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
MÉTODOS ABIERTOS Método de la Secante o de Newton Lagrange Método de Newton o de Newton Raphson Métodos de Puntos Fijos Año 2014
5 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE BISECCIÓN Paso Inicial Un intervalo [a0 ; b0] en el eje de las abscisas X en el cuál la función no lineal tenga al menos una raíz.
Esto es abscisas tales que f(a0).f(b0) < 0 Y(x)
f(b0) Y=f(x)
X
f(a0)
X= r a0
b0
Primera Aproximación Se debe encontrar la primera aproximación de la raíz; es decir, el primer elemento de la sucesión de soluciones aproximadas. Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
6 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE BISECCIÓN Primera aproximación Dado el Intervalo Inicial (a0; b0);
f(b0) se aproxima la raíz con El Punto Medio del Intervalo
r1
(a0 b0 ) 2
Y=f(x)
X
f(a0)
X= r a0
Se tienen dos intervalos
Intervalo (a0; r1)
r1
b0
Intervalo (r1; b0)
y se debe seleccionar uno de ellos que será el Intervalo 1 (a1; b1) en el cual está la raíz. ¿Cuál es?
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
7 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
Segunda aproximación
Dado el Intervalo 1 (a1; b1); f(b0) Y=f(x)
se aproxima la raíz con el Punto Medio del intervalo 1
r2
(a1 b1 ) 2
r1
X= r
X
f(a0) r1
a0 a1
r2
b0
b1
Intervalo (r2; b1) Se tienen dos intervalos: Intervalo (a1; r2) y se debe seleccionar uno de ellos que será el Intervalo 2 (a2; b2) en el cual está la raíz.
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
8 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
Tercera aproximación
Dado el Intervalo 2 (a2; b2); f(b0) Y=f(x)
se aproxima la raíz con el Punto Medio del intervalo 2
r3
(a2 b2 ) 2
X= r
X
f(a0) a0 a1 a2
r1 r2 r3
b0
b1
b2
Intervalo (r3; b2) Se tienen dos intervalos Intervalo (a2; r3) y se debe seleccionar uno de ellos que será el Intervalo 3 (a3; b3) en el cual está la raíz.
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
9 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
Cuarta aproximación
Dado el Intervalo 3 (a3; b3); se aproxima la raíz con el Punto Medio del intervalo 3
r4
f(b0)
Y=f(x)
(a3 b3 ) 2
X= r X
Se tienen dos intervalos Intervalo (a3; r4) Intervalo (r4; b3) y se debe seleccionar uno que será el Intervalo 4 (a4; b4) en el cual está la raíz. Dr. Ing. A.Mirasso
f(a0)
Raíces de Ecuaciones No Lineales
r1
a0 r2
a1 r3
a2 a3
Año 2014
b0
b1
b2 r4 b 3 10 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
Y(x)
k-ésima aproximación
Dado el Intervalo k (ak; bk) que cumple la condición de inicialización
Y=f(x) f(bk)
f(ak)*f( bk)<0 se calcula la nueva raíz aproximada como la abscisa media del Intervalo k
rk+1 = (ak + bk)/2 Se controla ! si
f(rk+1)
r
X
f(ak) ak
rk+1 es la solución buscada".
rk+1
bk
Se actualizan las variables si es que no se alcanzó la solución buscada, eligiendo de los dos intervalos definidos con ak; bk y rk+1. Si (f(ak)*f(rk+1)<0) entonces ak+1= ak bk+1 =rk+1
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Si (f(bk)*f(rk+1)<0) entonces ak+1= rk+1 bk+1 = bk
Año 2014
11 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
ALGORITMO DEL MÉTODO BISECCIÓN INICIALIZACIÓN
Con k=0, y el Intervalo (a0; b0) tal que f(a0)*f( b0)<0 No_Hay_Solu=Verdadero
Y(x) Y=f(x) f(bk)
HACER MIENTRAS (No_Hay_Solu) RECURRENCIA
rk+1 = (ak + bk)/2 CONTROL DE DETENCIÓN
f(rk+1)
r
X
f(ak)
Si ( f(rk+1 )=0) entonces No_Hay_Solu=Falso Fin SI
ak
rk+1
k =k+1 Si (f(ak)*f(rk+1)<0) entonces
ak+1= ak
bk =rk+1
Si
ak+1= rk+1
bk+1= bk
bk
ACTUALIZACIÓN
(f(bk)*f(rk+1)<0) entonces
FIN DEL HACER MIENTRAS Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
12 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE BISECCIÓN Ejemplo: Buscar una aproximación de
9
20 y=x^2-9
15 10 5 x
0 -6
-4
-2
0
2
4
6
-5 -10 -15
itera
a
b
f(a)
f(b)
r
f(r).
0
1
4
-8
7
2,5
-2,75
1 2 3
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
13 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
CONTROL DE DETENCIÓN Si ( f(rk+1 )=0) entonces No_Hay_Solu=Falso Fin SI Pero es !muy exigente" Es conveniente “suavizar”
Y(x)
Si (Valor absoluto de f(rk+1 )<İf) entonces No_Hay_Solu=Falso Fin SI Es conveniente !controlar el cambio de digitos" Si (Valor absoluto de [rk+1- rk ]<İ1) entonces No_Hay_Solu=Falso Fin SI
Raíces de Ecuaciones No Lineales
X= r X
Y(x)
Si (Valor absoluto de [(rk+1- rk)/ rk+1 ]<İr) entonces No_Hay_Solu=Falso Fin SI Dr. Ing. A.Mirasso
Y=f(x)
Año 2014
Y=f(x)
X
rk
rk+1
14 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE REGULA FALSI INICIALIZACIÓN Es igual que BISECCIÓN. Es decir, se necesita un intervalo donde al menos exista uan raíz
Y(x) Y=f(x) Recta L0
f(b0)
X b0
f(a0)
a0 Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
15 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
ALGORITMO DEL MÉTODO REGULA FALSI
Y(x)
INICIALIZACIÓN
Con k=0, y el Intervalo (a0; b0) tal que f(a0)*f( b0)<0 No_Hay_Solu=Verdadero
Y=f(x) Recta Lk
f(bk)
HACER MIENTRAS (No_Hay_Solu) RECURRENCIA rk 1
ak
f ( ak ) con mk
mk
§ f (ak ) f (bk ) · ¸¸ ¨¨ © (ak ) (bk ) ¹
CONTROL DE DETENCIÓN
bk
f(ak)
X= rk+1
Si ( f(rk+1 )=0) entonces No_Hay_Solu=Falso Fin SI
k =k+1 Si (f(ak)*f(rk+1)<0) entonces Si (f(bk)*f(rk+1)<0) entonces
X
ak
ACTUALIZACIÓN
ak+1= ak ak+1= rk+1
bk =rk+1 bk+1= bk
FIN DEL HACER MIENTRAS Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
16 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE REGULA FALSI Ejemplo: Buscar una aproximación de
9
20 y=x^2-9
15 10 5 x
0 -6
-4
-2
0
2
4
6
-5 -10 -15
itera
a
b
f(a)
f(b)
m
r
f(r).
0
1
4
-8
7
5
2,6
-2,24
1 2 3 4
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
17 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE LA SECANTE INICIALIZACIÓN Se necesita dos puntos cercanos a una raíz. Pueden ser dos raíces aproximadas con otro método
Y(x) Y=f(x)
f(rk-1)
f(rk)
X rk
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
rk-1
18 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
ALGORITMO DEL MÉTODO DE LA SECANTE INICIALIZACIÓN
Dos Puntos cercanos a la raíz buscada
Y(x) Y=f(x)
k=0, No_Hay_Solu=Verdadero
f(rk-1)
Recta Lk
HACER MIENTRAS (No_Hay_Solu) RECURRENCIA
rk 1
rk
f (rk ) con mk
§ f (rk ) f (rk 1 ) · ¸¸ ¨¨ © (rk ) (rk 1 ) ¹
mk
f(rk)
CONTROL DE DETENCIÓN Si ( f(rk+1 )=0) entonces No_Hay_Solu=Falso Fin SI
X rk+1
ACTUALIZACIÓN k =k+1
rk-1= rk rk= rk+1
rk
rk-1
Se retienen dos soluciones aproximadas
FIN DEL HACER MIENTRAS Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
19 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE LA SECANTE Ejemplo: Buscar una aproximación de
9
20 y=x^2-9
15 10 5 x
0 -6
-4
-2
0
2
4
6
-5 -10 -15
itera
a
b
f(a)
f(b)
m
r
f(r).
0
1
4
-8
7
5
2,6
-2,24
1 2 3 4
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
20 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE NEWTON RAPHSON INICIALIZACIÓN Se necesita un punto cercano a una raíz. Pueden ser una raíz aproximada con otro método, o simplemente una propuesta
Y(x) Y=f(x)
f(rk)
rk
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
X
21 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
ALGORITMO DEL MÉTODO NEWTON RAPHSON INICIALIZACIÓN
Un Punto cercano a la raíz buscada k=0, No_Hay_Solu=Verdadero
Y(x)
HACER MIENTRAS (No_Hay_Solu)
Y=f(x)
Tk
f(rk)
RECURRENCIA
rk 1
rk
f (rk ) con mk
df dx
mk
rk
CONTROL DE DETENCIÓN Si ( f(rk+1 )=0) entonces No_Hay_Solu=Falso Fin SI
rk+1
ACTUALIZACIÓN k =k+1
rk= rk+1
X
rk
Se retienen una solución aproximada
FIN DEL HACER MIENTRAS
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
22 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE NEWTON RAPHSON Ejemplo: Buscar una aproximación de
9
20 y=x^2-9
15 10 5 x
0 -6
-4
-2
0
2
4
6
-5 -10 -15
itera
rk
f(rk)
m
0
1
-8
2
rk+1
f(rk+1).
1 2 3 4 5
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
23 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE PUNTO FIJO
y=F(x)
Dada una función no lineal
F(x)
y = F(x) se busca el valor de abscisa xs tal que
C
F ( xs ) C Con C una constante arbitraria y conocida.
xs
Es posible escribir la ecuación no lineal en la forma:
\ ( xs )
F ( xs ) C
y \ ( x)
F ( x) C
<(x)
0
y=<(x) y=F(x)
que es la búsqueda de la raíz de la función
El valor de la abscisa xs es INVARIANTE a que la ecuación no lineal se multiplique por D un ESCALAR NO NULO
X
Y=D . <(x)
F(x)
X
C
D \ ( x s ) D ( F ( x s ) C ) 0
xs
X
Es posible sumar en ambos miembros el valor de la abscisa xs Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
24 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
D \ ( x s ) D ( F ( x s ) C ) 0
Si a la ecuación no lineal
<(x)
xs D \ ( xs )
y=x
y=<(x)
Se le suma xs en ambos miembros resulta
xs
que se puede escribir en la forma:
xs donde
g ( x)
y=F(x)-C
y=x + D . <(x)
F(x)
Y=D . <(x)
g ( xs )
x D ( F ( x) C )
X
C
x D \ ( x)
xs
X
La igualdad x g (x) se puede interpretar como la intersección de las siguientes funciones Recta que bisecta el primer cuadrante (y=x) con la curva y=g(x). El punto solución es tal que tiene abscisa y ordenadas de igual valor. Es el Punto Fijo de a curva g(x)
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
Año 2014
25 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE PUNTO FIJO
Inicialización Se necesita un punto cercano a una raíz. (Igual a Newton Raphson) <(x)
Recurrencia La aproximación de la raíz se obtiene mediante,
xk 1
g ( xk )
Y=F(x) - C
y=x+D . <(x)
g(xk)
Control de Detención Igual que métodos anteriores. Alternativamente se debe controlar si se cumple que xk 1 g ( xk 1 ) d H f
g(xk+1) y=g(x)
g(xs) F(x)
X xs
C
Actualización de Variables Se deben retener la última aproximación obtenida. Dr. Ing. A.Mirasso
y=x
Raíces de Ecuaciones No Lineales
xk+1
xr X
Año 2014
26 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE PUNTO FIJO
<(x)
<(x)
y=x
g(x2) g(x1)
Y=F(x) - C
g(x0) g(x1) g(xs)
y=g(x)
y=x y=F(x) - C
g(x0) g(xs) F(x)
F(x)
C
y=g(x)
xs
x2
x1
x0
X xs
C
x1
x2 x3
x0 X
X
CONVERGE Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
X
DIVERGE Año 2014
27 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE PUNTO FIJO Condición de Convergencia del Método de Punto Fijo
La solución de la igualdad
xs
x
g (x )
es xs tal que es el punto fijo de g(x); es decir, <(x)
g ( xs )
y=x Y=F(x) - C
La fórmula de recurrencia es
xk 1
g ( xk )
xk 1 xs
g ( xk ) g ( x s )
Al restar miembro a miembro, se tiene
g(x0) g(x1) g(xs)
y=g(x)
F(x)
Al aplicar el Teorema del Valor Medio, resulta ( xk 1 xs )
dg ( x ) ( xk x s ) dx x [
xs
C
con [ una abscisa entre xk y xs..
H k 1
dg ( x ) dx
Dr. Ing. A.Mirasso
x [
Hk
CONVERGE si Raíces de Ecuaciones No Lineales
x2
x1
x0
X
X dg ( x) dx x Año 2014
[
1 28 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE PUNTO FIJO
9 , es buscar que x2=9.
Ejemplo: Buscar una aproximación de y=x^2-9
9
g(x)
7
y=x
rk
rk+1=g(rk)
f(rk+1).
0
1
1,8
-5,76
1
5 3
2
1 -1
itera
0
2
4
6
8
10
x 12
3
-3
4
-5
5
-7 -9
6
-11
7
g ( x) dg ( x) dx Dr. Ing. A.Mirasso
x (1 / 10) ( x 2 9)
8
1 (2 / 10) x
Raíces de Ecuaciones No Lineales
Año 2014
29 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO ITERATIVOS PARA OBTENER RAICES DE ECUACIONES NO LINEALES
Método Bisección Regula Falsi Secante
Newton Raphson
Resumen de los Métodos Inicialización Recurrencia Intervalo con rk+1 = (ak + bk)/2 Raíz Intervalo con § f (ak ) f (bk ) · f (ak ) ¸¸ con mk ¨¨ rk 1 ak Raíz mk © (ak ) (bk ) ¹ 2 Puntos Cercanos 1 Punto Cercano
Punto Fijo 1 Punto Cercano
Dr. Ing. A.Mirasso
rk 1
rk 1
xk 1
f (rk ) con rk mk rk
f (rk ) con mk
mk
mk
§ f (rk ) f (rk 1 ) · Retiene 2 ¸¸ ¨¨ © (rk ) (rk 1 ) ¹ Puntos df dx
rk
Cercanos Retiene 1 Punto Cercano Retiene 1 Punto Cercano
g ( xk )
Raíces de Ecuaciones No Lineales
Actualización Elije Intervalo con Raíz Elije Intervalo con Raíz
Año 2014
30 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO ITERATIVOS PARA OBTENER RAICES DE ECUACIONES NO LINEALES
CONTROL DE DETENCIÓN Los siguientes controles son complementarios
Y(x)
Si (Valor absoluto de f(rk+1 )<İf) entonces No_Hay_Solu=Falso Fin SI
Y=f(x) X= r X
Si (Valor absoluto de [rk+1- rk ]<İ1) entonces No_Hay_Solu=Falso Fin SI
Y(x)
Y=f(x)
Si (Valor absoluto de [(rk+1- rk)/ rk+1 ]<İr) entonces No_Hay_Solu=Falso Fin SI Si (k > MaxIter) entonces No_Hay_Solu=Falso Fin SI
Dr. Ing. A.Mirasso
Raíces de Ecuaciones No Lineales
X
rk
Año 2014
rk+1
31 de 31
Facultad de Ingeniería, Universidad Nacional de Cuyo
Introducción. Repaso de métodos iterativos para obtener raíces de ecuaciones no lineales Sistemas de Ecuaciones Lineales (No Homogéneos) Método de Jacobi. Algoritmo, Convergencia Método de Gauss Seidel. Algoritmo Autovalores y Autovectores Propiedades Método de la Potencia. Algoritmo, Convergencia Método de la Potencia Inversa El cociente de Rayleigh Aplicaciones
Facultad de Ingeniería, Universidad Nacional de Cuyo
Los métodos iterativos generan una sucesión de soluciones aproximadas, que debe converger a la solución exacta del problema que se pretende resolver. Cada iteración genera una solución aproximada, y su diferencia respecto de la solución exacta debe ser menor que en la iteración anterior para que exista convergencia Tienen en general los siguientes componentes
Inicialización o Acercamiento Recurrencia Control de Detención Actualización y deben cumplir algún CRITERIOS DE CONVERGENCIA Para garantizar que se encontrará la solución del problema
Facultad de Ingeniería, Universidad Nacional de Cuyo
Se busca el valor Síntesis:
xk +1 = r
Punto Fijo
Inicialización Recurrencia
f (r ) = 0
tal que
Método Iterativo Newton Raphson
x0
Se propone una solución inicial
xk +1 = g ( xk )
rk +1 = rk −
f ( xk +1 ) ≤ ε f
Control de Detención Actualización CRITERIOS DE CONVERGENCIA
xk +1 − xk ≤ ε a
mk =
df dx
rk
k ≤ k max
o o
f (rk ) con mk
xk +1 − xk ≤ εr xk +1
Retiene la última solución aproximada dg ( x ) <1 dx x=ξ
dg ( x) dx
< 1 con g ( x ) = x − x =ξ
f ( x) df dx
Facultad de Ingeniería, Universidad Nacional de Cuyo
La implementación se puede sintetizar de la siguiente forma: x0 ε a ε r ε f Se eligen los valores de Inicialización Se define ban =true while(ban) f (rk ) r r = − Recurrencia k +1 k m con k
xk +1 = g ( xk ) k=k+1 if
f ( xk +1 ) ≤ ε f
(
)
ban =false
Control de Detención
end if
(
xk +1 − xk ≤ ε a )
ban =false end if ( k ≤ k max )
ban =false end
xk = xk +1 end
Actualización de Variables
mk =
df dx
rk
Facultad de Ingeniería, Universidad Nacional de Cuyo
Dada una Matriz A ε RNxN, un vector b ε RN, el sistema de ecuaciones lineales asociado es
A⋅ x = b b1
a11
a12
a13
....
a1N
x1
a21 a31
a 22 a32
a 23 a33
.... ....
a2 N a3 N
b2 x2 x3 = b3
.... aN1
.... aN 2
.... aN 3
.... ....
.... a NN
... xN
... bN
A = D+B
Es posible escribir a11 0 0 ....
0 a 22 0 ....
0 0 a33 ....
.... .... .... ....
0 0 0 ....
0
0
0
....
a NN
=D
0 a 21 B = a31 ....
a12 0 a32 ....
a13 a 23 0 ....
.... .... .... 0
a1N a2N a3 N ....
a N1
aN 2
aN3
....
0
El sistema se puede escribir
D⋅ x + B⋅ x = b
Facultad de Ingeniería, Universidad Nacional de Cuyo
Dada una Matriz A ε RNxN, un vector b ε RN, el sistema de ecuaciones lineales asociado es
A⋅ x = b El sistema se puede escribir
D⋅ x + B⋅ x = b x = −D−1 ⋅ B ⋅ x + D−1 ⋅ b
x = T⋅ x + c T = −D−1 ⋅ B 0 − a 21 / a 22 T = − a31 / a33 .... − a N 1 / a NN
− a12 / a11 0 − a32 / a33 .... − a N 2 / a NN
− a13 / a11 − a 23 / a 22 0 .... − a N 3 / a NN
c = D−1 ⋅ b .... .... .... 0 ....
− a1N / a11 − a 2 N / a 22 − a3 N / a33 .... 0
Se puede iterar con
x
( k +1)
= T⋅ x
(k )
+c
Hasta encontrar que el ERROR es tan pequeño como se quiera
b1 / a11 b2 / a 22 c = b3 / a33 ... bN / a NN
Facultad de Ingeniería, Universidad Nacional de Cuyo
EJEMPLO Se busca la solución del siguiente sistema de ecuaciones lineales
− 25
x1 10 x2 − 25 50 − 25 20 0 = 20 0 − 25 50 − 25 x3 10 0 0 − 25 50 x4 50
0
0
A⋅ x = b
La matriz de coeficientes se puede escribir
=
A 50
− 25
0
0
+
D 50
0 − 25 50 − 25 0 = 0 0 − 25 50 − 25 0 0 − 25 50 0
0
0
0
B 0
− 25
0
0
50 0 0 − 25 0 − 25 0 + 0 50 0 0 0 − 25 − 25 0 0 50 0 0 − 25 0
Facultad de Ingeniería, Universidad Nacional de Cuyo
x = −D−1 ⋅ B ⋅ x + D−1 ⋅ b
x = T⋅ x + c D −1
c=
c=
⋅
b
50 0 0 50
0 0
0 0
50 0 0 50
0 0
−1
10 20
⋅
20 10
− D −1
T=
0 T=−
0 0
0 0
0
0
0
=
0 0 0 0 0 0 0 0 0 0 0 0
⋅
⋅
10 20
=
20 10
B
0
0
0
− 25
0
0
0
0
− 25
0
− 25
0
0
− 25
0
− 25
0
0
− 25
0
0 0
⋅
=
Facultad de Ingeniería, Universidad Nacional de Cuyo
A⋅ x = b 50
− 25
0
0
x1
10
− 25
50
− 25
0
x2
20
0
− 25
50
− 25
0
0
− 25
50
x3 x4
=
20 10
x = −D −1 ⋅ B ⋅ x + D −1 ⋅ b x= T ⋅x + c x1 0 25 / 50 0 0 10 / 50 x1 25 / 50 0 25 / 50 0 20 / 50 ⋅ + x= x3 0 25 / 50 0 25 / 50 20 / 50 x4 0 0 25 / 50 0 10 / 50
Facultad de Ingeniería, Universidad Nacional de Cuyo
Proceso Iterativo de JACOBI
x
x
( k +1)
( k +1)
= T⋅ x
(k )
+c
x1 0 25 / 50 0 0 x2 25 / 50 0 25 / 50 0 = ⋅ x3 0 25 / 50 0 25 / 50 x4 0 0 25 / 50 0
(k )
10 / 50 20 / 50 + 20 / 50 10 / 50
(0)
Para un vector inicial x arbitrario Iter. x1 x2 0 1,5 2,5 (25/50)*(2,5)+10/50 (25/50)*(1,5)+(25/50)*(2,5)+20/50 1 1 1,45 2,4 2 3 4 Se puede paralelizar ¡!
x3
x4
2,5
1,5
(25/50)*(2,5)+(25/50)*(1,5)+20/50 (25/50)*(2,5)+10/50
2,4
1,45
Facultad de Ingeniería, Universidad Nacional de Cuyo
Proceso Iterativo de JACOBI
x
x
( K +1)
( k +1)
= T⋅ x
(k )
+c
x1 0 25 / 50 0 0 x1 25 / 50 0 25 / 50 0 ⋅ = 0 25 / 50 0 25 / 50 x3 x4 0 0 25 / 50 0
(K )
10 / 50 20 / 50 + 20 / 50 10 / 50
Iter.
x1
x2
x3
x4
dif_x1 dif_x2 dif_x3 dif_x4 Norma dif
0 1
1,5
2,5
2,5
1,5
1,45
2,4
2,4 2,325
1,45
------1,45 -1,5
2 3
1,3625
Norma de un vector: Cuadrática o Infinito.
------2,4 -2,5
--------0,1
--------0,05
-------0,1 0,075 0,0625
Facultad de Ingeniería, Universidad Nacional de Cuyo
A partir del método iterativo de Jacobi, cuyo fórmula de recurrencia está dada por
x
( k +1)
= T⋅ x
( k +1)
(k )
+c
( k +1)
(k )
x = Tl ⋅ x + Ts ⋅ x Se puede iterar con Hasta encontrar que el ERROR es tan pequeño como se quiera. Siendo 0
0
0
....
0
x1( k +1)
− a 21 / a 22
0
0
....
0
x 2( k +1)
Tl = − a31 / a33 ....
− a32 / a33 ....
0 ....
.... 0
0 ....
− a N 1 / a NN
− a N 2 / a NN
− a N 3 / a NN
....
0
x N( k +1)
− a12 / a11 0 0 .... 0
− a13 / a11 − a 23 / a 22 0 .... 0
− a1N / a11 − a 2 N / a 22 − a 3 N / a33 .... 0
x1( k ) x 2( k )
Ts =
0 0 0 .... 0
.... .... .... 0 ....
x
( k +1)
+c
= x3( k +1) ...
x
(k )
= x3( k ) ... x N( k )
Se debe destacar que: • para calcular x1, participa todo el vector x de la iteración anterior. • Para calcular x2, participa x1 de la iteración actual (recién calculado)y todas las demás componentes del vector x de la iteración anterior. • Para calcular x3, participa x1 y x2 de la iteración actual (recién calculadas)y todas las demás componentes del vector x de la iteración anterior. • Para calcular x4, participa x1 ,x2 y x3 de la iteración actual (recién calculadas)y todas las demás componentes del vector x de la iteración anterior. • Así se sigue hasta calcular xN con todas las componentes de la iteración actual del vector x (recién calculadas), desde la 1 hasta la N-1.
Facultad de Ingeniería, Universidad Nacional de Cuyo
x
Se itera con
= Tl ⋅ x
( k +1)
+ Ts ⋅ x
(k )
+c
0
0
0
....
0
x1( k +1)
− a 21 / a 22
0
0
....
0
x 2( k +1)
− a 32 / a33
0
....
0
....
....
....
0
....
...
− a N 1 / a NN
− a N 2 / a NN
− a N 3 / a NN
....
0
x N( k +1)
Tl = − a 31 / a33
Ts =
( k +1)
x
0
− a12 / a11
− a13 / a11
....
− a1N / a11
0
0
− a 23 / a 22
....
− a 2 N / a 22
0
0
0
....
− a3 N / a33
.... 0
.... 0
.... 0
0 ....
.... 0
( k +1)
= x3( k +1)
x1( k ) x 2( k ) x
(k )
= x3( k ) ... x N( k )
Se debe destacar que: • para calcular x1, participa todo el vector x de la iteración anterior. • Para calcular x2, participa x1 de la iteración actual (recién calculado) y todas las demás componentes del vector x de la iteración anterior. • Así se sigue hasta calcular xN con todas las componentes de la iteración actual del vector x (recién calculadas), desde la 1 hasta la N-1.
Facultad de Ingeniería, Universidad Nacional de Cuyo
Síntesis Punto Fijo Se busca
Inicialización Recurrencia
Control de Detención Actualización CRITERIOS DE CONVERGENCIA
Sistema de Ecuaciones Lineales
xk +1 = r tal que f (r ) = 0
x ε R N tal que
x0
Se propone una solución inicial
xk +1 = g ( xk )
x
f ( xk +1 ) ≤ ε f
( k +1)
o
= T⋅ x
f = A⋅x −b = 0
(k )
+c
x
xk +1 − xk ≤ ε a
( k +1)
= Tl ⋅ x
( k +1)
Ts ⋅ x o
(k )
+
+c
k ≤ k max
Retiene la última solución aproximada dg ( x ) <1 dx x=ξ
En lugar de Valor absoluto de una variable real
Mayor de los Valores Matriz A Absolutos de T estrictamente debe ser menor a 1 diagonal dominante ρ (T) < 1 se usa Norma de un Vector de componentes reales
Facultad de Ingeniería, Universidad Nacional de Cuyo
Dada una matriz A, se le puede asociar una Transformación Lineal, la que tiene Direcciones Invariantes, soluciones de
(A − λ ⋅ I) x = 0 Ejemplos Simetría respecto de una RECTA
Proyección en EJE X en dirección de una RECTA
y
y Tu
u u x Dirección Invariante 1 Autovalor Dirección Invariante 2 Autovalor
Tu x
Facultad de Ingeniería, Universidad Nacional de Cuyo
Dada una matriz A, se buscan las Direcciones Invariantes, soluciones de
(A − λ ⋅ I) x = 0 o bien A ⋅ x = λ ⋅ x Esto es equivalente a encontrar un vector
y
y = A⋅ x
que se puede obtener como
y = λ⋅x
y resulta igual a lo que significa que los vectores El METODO DE LA POTENCIA es un algoritmo iterativo que propone un Obtiene un nuevo vector como Sólo se detiene si se cumple que EJEMPLO: Dada la matriz A -10 -4
y0 2 2 iterac alfa(1) alfa(2)
y1 -12 -8 1
y0
x
e
y
son PARALELOS
arbitrario y
yk +1 = A ⋅ yk yk +1
es paralelo a
yk
4 0
y2 88 48 2
y3 -688 -352 3
y4
y5
y6
4
5
6
Facultad de Ingeniería, Universidad Nacional de Cuyo
El METODO DE LA POTENCIA es un algoritmo iterativo para resolver
y0
que propone un
arbitrario y se escala para obtener el versor
yk +1
Sólo se detiene si se cumple que EJEMPLO: Dada la matriz A
y0 2 2 norma
es paralelo a
xk
4 0
y1 -6 -4
y2
y3
y4
y5
y6
x0 x1 1 1 1 0,66666667
x2
x3
x4
x5
x6
2
3
4
5
6
2
-6
iterac alfa(1) alfa(2) num den rho
x0 , con el que se inicia el proceso iterativo
yk +1 = A ⋅ xk
Obtener un nuevo vector como
-10 -4
A⋅ x = λ ⋅ x
1
Facultad de Ingeniería, Universidad Nacional de Cuyo
Se busca
x ε R N tal que
Inicialización Recurrencia Control de Detención
Autovalores
Sistema de Ecuaciones Lineales
f = A⋅ x −λ ⋅ x = 0
f = A⋅x −b = 0
Se propone una solución inicial
y k +1 = A ⋅ xk
(
x k +1 = y k +1 ⋅ 1 / y k +1
f ( xk +1 ) ≤ ε f
Actualización CRITERIOS DE CONVERGENCIA
o
x0 x
)
( k +1)
xk +1 − xk ≤ ε a
= T⋅ x o
(k )
+c
k ≤ k max
Retiene la última solución aproximada La matriz A debe ser Diagonalizable. Sus autovectores deben ser linealmente independientes
En lugar de Valor absoluto de una variable real
Mayor de los Valores Absolutos de T debe ser menor a 1 ρ (T) < 1 Matriz A estrictamente diagonal dominante
se usa Norma de un Vector de componentes reales
Facultad de Ingeniería, Universidad Nacional de Cuyo
Sea el sistema EDO
y (t ) = A ⋅ y (t ) admite como solución y (t ) autovectores
con
= v ⋅ e λ ⋅t siendo
v, λ
v − 10 4 v1 =λ⋅ 1 v2 − 4 0 v2
La solución general será la combinación lineal λ 1⋅t λ 2⋅t 1 2
y (t ) = v ⋅ e
5 y ( 0) = 3
+ v ⋅e
,
− 10 4
A=
−4
0
la solución del siguiente sistema de autovalores y 6 y1(t)
5
y2(t)
4 3 2 1 0 0
0,5
1
1,5
1 1 −2⋅t 14 1 −8⋅t = ⋅ + t e ⋅ e y ( ) Es decir, 3 2 3 0.5
y2
que cumple con las condiciones iniciales. 3,5 3 2,5 2 1,5 1 0,5
y1
0 0
1
2
3
4
5
6
Facultad de Ingeniería, Universidad Nacional de Cuyo
Dada una matriz A ε RNxN, tal que Sus autovectores v1, v2, …vN, y sus autovalores λ1, λ2, λ3,…. λN, verifican que
A . vk= λk . vk
con k=1 a N
Es diagonalizable: es decir tiene N autovectores linealmente independientes que forman base de RN. Es posible ordenar los autovalores de mayor a menor considerando sus valores absolutos en la forma:
λ1 > λ2 > λ3 > ........... λ N Entonces el algoritmo del método de la potencia converge
Facultad de Ingeniería, Universidad Nacional de Cuyo
Se busca u(x,t) solución de
u(x,t) U1(t)
2 ∂ 2 u ( x, t ) 2 ∂ u ( x, t ) − 12 +x = 0 en Ω = {x ε R : 0 ≤ x ≤ 1} ∂x 2 ∂t 2 u (0, t ) = 0 u (1, t ) = 0
U0(t)
0,25
U2(t)
0,5
U3(t) 0,75
U4(t) X=1
Se plantea encontrar una solución aproximada en forma de función discreta en la variable x, t aunque continua en la variable t. Se pretende encontrar las funciones Uk(t)=u(Xk,t) con k=0,N, en N+1 puntos elegidos del dominio x, identificados con su abscisa Xk.
M ⋅ u (t ) + K ⋅ u(t ) = 0 con
2
−1
0
12 K= −1 2 −1 ; M = 0,25 2 0 −1 2
0,25 2 0 0
0 0,50 0
0 2
Que admite una solución del tipo
u (t ) = x ⋅ e I ω t Que resulta en
⋅ (K − ω 2 M ) ⋅ x = 0
0 0,75 2
U 1 (t ) u(t ) = U 2 (t ) U 3 (t )
X
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN Y APROXIMACIÓN DE FUNCIONES
Planteo de Problema Interpolación de funciones Método Directo Método con Polinomios de Lagrange Método con Polinomios de Newton Aproximación de funciones Método de Mínimos Cuadrados Ejemplos Resumen
Dr. Ing. A.Mirasso
Interpolación y Aproximación
Año 2014
1 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN Y APROXIMACIÓN DE FUNCIONES DISCRETAS
Dada una función discreta, y=f(x) de R o R definida mediante (n+1) puntos (xi; yi=f(xi)) con i=0,n. con k=0,m. Se propone Pm(x)= 6 ak Ik(x) Donde ak son las incógnitas; y {Ik(x)} es una Base elegida Se define: Residuo ri = f(xi) - Pm(xi) con k=0,m; i=0,n r0 ½ °r ° °° 1 °° ®r2 ¾ °... ° ° ° °¯rn °¿
ri =
yi
y 0 ½ § I 0 ( x0 ) ° y ° ¨ I (x ) °° 1 °° ¨ 0 1 ¨ ® y 2 ¾ ¨ I0 ( x2 ) ° ... ° ¨ .... ° ° ¨ °¯ y n °¿ © I 0 ( x n )
-
6 ak Ik(xi)
I1 ( x0 ) I1 ( x1 ) I1 ( x 2 )
I 2 ( x0 ) I 2 ( x1 ) I2 ( x2 )
....
....
I1 ( x n )
I2 ( xn )
.... .... .... .... ....
x
I m ( x0 ) · a 0 ½ ¸ I m ( x1 ) ¸°° a1 °° ° ° I m ( x 2 ) ¸® a 2 ¾
x0
¸ .... ¸° ... ° ° ° I m ( x n ) ¸¹°¯a m °¿
INTERPOLACIÓN Forma Fuerte: ri=0, para todo xi f(x)
y
y
xi
APROXIMACIÓN Forma Débil: los ri !son nulos en promedio f(x)
y
ri
Dr. Ing. A.Mirasso
xn Interpolación y Aproximación
r1
r0
x xi
y )a
r
Pm(x)
x0
xn
r2
x x0
Año 2014
rn Pm(x)
xi
xn 2 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN: MÉTODO DIRECTO
Dada una función discreta, y=f(x) de R o R definida mediante (n+1) puntos (xi; yi=f(xi)) con i=0,n.
Se adopta como Base = {1, x ,x2, x3,......,xn}
Pn(x)= 6 (ak x ) k
Se propone
Se determinan los ak tales que el Residuo sea nulo, es decir
r0 ½ °r ° °° 1 °° ®r2 ¾ °...° ° ° °¯rn °¿
r
y0 ½ § I0 ( x0 ) ° y ° ¨ I (x ) °° 1 °° ¨ 0 1 ¨ ® y2 ¾ ¨ I0 ( x2 ) ° ... ° ¨ .... ° ° ¨ °¯ yn °¿ © I0 ( xn )
y )a
I1 ( x0 ) I1 ( x1 ) I1 ( x2 )
I2 ( x0 ) I2 ( x1 ) I2 ( x2 )
.... .... ....
.... I1 ( xn )
.... I2 ( xn )
.... ....
§ 1 x0 x ¨ ¨ 1 x1 x ¨1 x x 2 ¨ ¨ .... .... .... ¨ 2 © 1 xn xn
¸ .... ¸° ... ° ° ° Im ( xn ) ¸¹°¯am °¿
f(x)
0 ½ °0 ° °° °° ®0 ¾ °...° ° ° °¯ 0 °¿
Pn(x)
x x0
xi
xn
.... x ·a0 ½ ¸° ° .... x ¸° a1 ° ° ° .... x ¸¸®a 2 ¾ .... .... ¸° ... ° ° ° n¸ .... x n ¹°¯a n °¿ n 0 n 1 n 2
I1(xi)
I2(xi)
y
I0(xi)
y0 ½ °y ° °° 1 °° ® y2 ¾ ° ... ° ° ° °¯ y n °¿
x x0
para que sea Solución Única; m=n y
Resulta el polinomio interpolante Dr. Ing. A.Mirasso
con k=0,n
0
Que resulta en el sistema de ecuaciones lineales 2 0 2 1 2 2
Im ( x0 ) · a0 ½ ¸ Im ( x1 ) ¸°° a1 °° ° ° Im ( x2 ) ¸® a2 ¾
y
Interpolación y Aproximación
xi
xn
todos x distintos!!!!!
Pn(x)= 6 (ak xk ) con k=0,n Año 2014
3 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN: MÉTODO DIRECTO
EJEMPLO
Dada una función discreta, y=f(x) de R o R definida mediante (2) puntos (xi; yi=f(xi)) con i=0,1.
x 1,5 3
y=f(x) 3 4
y P1(x) f(x)
Se adopta como Base = {1, x} con lo que la matriz ) y el sistema a resolver resultan:
§ ¨¨ ©
·a0 ½ ¸¸® ¾ ¹¯ a1 ¿
de donde se tiene que a0
½ ® ¾ ¯ ¿
x0 y
x x1
El Polinomio Interpolante usando la Base de Polinomios elementales {1,
P1 ( x)
Dr. Ing. A.Mirasso
1
I1(xi)
I0(xi)
x0
a1
x1
x} es
x
Interpolación y Aproximación
Año 2014
4 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN: MÉTODO DIRECTO
EJEMPLO
Dada una función discreta, y=f(x) de R o R definida mediante (3) puntos (xi; yi=f(xi)) con i=0,2. y
x y=f(x) 1,5 3 3 4 4,5 3,5 Se adopta como Base = {1, x, x2 }
f(x)
Pn(x)
con lo que la matriz ) y el sistema a resolver resultan:
§ ¨ ¨ ¨ ©
· a 0 ½ ¸° ° ¸® a1 ¾ ¸ °a ° ¹¯ 2 ¿
° ® ° ¯
½ ° ¾ ° ¿
x0
x1
y
I2(xi) I0(xi) x x0
x1
x2
a1 a2 de donde se tiene que a0 El Polinomio Interpolante usando la Base de Polinomios elementales {1,
P2 ( x) Dr. Ing. A.Mirasso
1
x
Interpolación y Aproximación
x2
Año 2014
I1(xi)
x2
x, x2} es 5 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN: MÉTODO DE POLINOMIOS DE LAGRANGE
Dada una función discreta, y=f(x) de R o R definida mediante (n+1) puntos (xi; yi=f(xi)) con i=0,n.
Se adopta como Base = {l0(x), l 1(x), l 2(x),...... ln(x)}, Los polinomios de Lagrange l i(x) se definen como, li (x)
(x x 1 )(x x 2 )(x x 3 )...........(x x n ) li (x) , (x i x 1 )(x i x 2 )(x i x 3 )...........(x i x n )
n
(x x k ) (x i x k ) . 0
que el Residuo sea nulo,
De donde resulta
y )a
0
§ 1 0 0 .... 0 ·a0 ½ ¸° ° ¨ 0 1 0 .... 0 ¸° a1 ° ¨ ¨ 0 0 1 .... 0 ¸°a ° ¸® 2 ¾ ¨ ¨ .... .... .... .... .... ¸° ... ° ¨ 0 0 0 .... 1 ¸°°a °° ¹¯ n ¿ ©
y0 ½ °y ° °° 1 °° ® y2 ¾ ° ... ° ° ° °¯ yn °¿
y resulta el polinomio interpolante Interpolación y Aproximación
Pn(x)
l0(x)
NO HAY QUE RESOLVER SISTEMA DE ECUACIONES!!!!
Dr. Ing. A.Mirasso
f(x)
k k zi
Son tales que para los i, k abscisas datos li(xi)=1 y li(xk)=0 con kzi Se determinan los ak imponiendo
r
y
Año 2014
lk(x)
x0
x0
I0(xi) x
xi
xn
Ik(xi) x
x0
xn
xi
Pn(x)= 6 (yk lk(x)) con k=0,n 6 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN: MÉTODO DE POLINOMIOS DE LAGRANGE
EJEMPLO
Dada una función discreta, y=f(x) de R o R definida mediante (2) puntos (xi; yi=f(xi)) con i=0,1.
x 1,5 3
y
y=f(x) 3 4
f(x)
x
l 0 ( x) y
Base = {l0(x), l1(x)} con
§ ¨¨ ©
½ ® ¾ ¯ ¿
de donde se tiene que a0
x0
x1
l0(x)
1
l1 ( x)
y así la matriz ) y el sistema a resolver resultan:
·a0 ½ ¸¸® ¾ ¹¯ a1 ¿
P1(x)
x y
x0
x1 l1(x)
1
x
x0
a1
x1
El Polinomio Interpolante usando la Base de Polinomios de Lagrange es
P1 ( x) Dr. Ing. A.Mirasso
Interpolación y Aproximación
Año 2014
7 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN: MÉTODO DE POLINOMIOS DE LAGRANGE
EJEMPLO
Dada una función discreta, y=f(x) de R o R definida mediante (3) puntos (xi; yi=f(xi)) con i=0,2.
x 1,5 3 4,5
y=f(x) 3 4 3,5
l 0 ( x) l1 ( x) l 2 ( x)
Con la Base = {l0(x), l1(x), l2(x)} resulta § ¨ ¨ ¨ ©
· a0 ½ ¸° ° ¸® a1 ¾ ¸ °a ° ¹¯ 2 ¿
° ® ° ¯
½ ° ¾ ° de donde se tiene que a0 ¿
a1
a2
El Polinomio Interpolante usando la Base de Polinomios de Lagrange { l0(x), l1(x), l2(x)} es
P2 ( x) AL AGREGAR UN PUNTO SE DEBE HACER TODO DE NUEVO Dr. Ing. A.Mirasso
Interpolación y Aproximación
Año 2014
8 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN: MÉTODO DE POLINOMIOS DE NEWTON
Dada una función discreta, y=f(x) de R o R definida mediante (n+1) puntos (xi; yi=f(xi)) con i=0,n. Se toma como base a los llamados polinomios de Newton. Estos tienen como particularidad que se basan en los polinomios bases anteriores.
y
n0(x) = 1 n1(x) = n0(x) (x-x0)
n0(x) = 1 n1(x) = 1(x-x0)
n2(x) = n1(x) (x-x1)
n2(x) = 1(x-x0) (x-x1)
n3(x) = n2(x)
n3(x) = 1(x-x0)(x-x1)(x-x2)
(x-x2)
f(x) Pn(x)
nk(x) = nk-1(x)·(x-xk-1), para todo kt1 Base = {n0(x),
Con la
n1(x), n2(x),...... nn(x)},
Se determinan los ak tales que el Residuo sea nulo, y resulta
§1 ¨ ¨1 ¨1 ¨ ¨ .. ¨1 ©
0 x1 x0
0 0
x 2 x0
( x 2 x0 )( x 2 x1 )
...........
...........................
x n x0
De donde:
( x n x0 )( x n x1 )
a0= y0,
Resulta el polinomio interpolante Dr. Ing. A.Mirasso
x
r
y )a
0
·a 0 ½ ¸° ° ¸° a1 ° ¸ °a ° ........ 0 ¸® 2 ¾ ....... 0 ¸° ... ° ° ° ........ ( x n x0 )( x n x1 ).....( x n x n 1 ) ¸¹°¯a n °¿
......... .........
0 0
a1
y1 y 0 , x1 x 0
a2
Pn(x)= 6 (ak nk(x))
Interpolación y Aproximación
n2(x)
n1(x)
y y0 ½ °y ° °° 1 °° ® y2 ¾ ° ... ° ° ° °¯ y n °¿
xn
n0(x) x x0
x1
x2
y 2 y1 y1 y 0 x 2 x1 x1 x 0 x2 x0
Año 2014
con k=0,1,2,!n 9 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN: MÉTODO DE POLINOMIOS DE NEWTON
EJEMPLO
Dada una función discreta, y=f(x) de R o R definida mediante (2) puntos (xi; yi=f(xi)) con i=0,1.
x
y=f(x)
1,5 3
3 4
'
'
y
f(x)
Base = {n0(x), n1(x)} con
n0 ( x) 1
x
y así la matriz ) y el sistema a resolver resultan:
§ ¨¨ ©
P1(x)
n1 ( x)
· a 0 ½ ½ ¸¸® ¾ ® ¾ ¹¯ a1 ¿ ¯ ¿
de donde se tiene que
a0
y
x0
x1
n0(x) n1(x)
1
x
x0
x1
a1
El Polinomio Interpolante usando la Base de Polinomios de Newton es
P1 ( x) Dr. Ing. A.Mirasso
Interpolación y Aproximación
Año 2014
10 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTERPOLACIÓN: MÉTODO DE POLINOMIOS DE NEWTON
EJEMPLO
Dada una función discreta, y=f(x) de R o R definida mediante (3) puntos (xi; yi=f(xi)) con i=0,2.
x 1,5 3 4,5
y=f(x)
'
'
y f(x)
3 4 3,5
Pn(x)
x
Base = {n0(x), n1(x), n2(x)} n0 (x) con
n2(x)
1
xn n1(x)
y
n0(x)
n1 ( x )
n2 (x) la ) y el sistema a resolver resultan:
§ ¨ ¨ ¨ ©
Dr. Ing. A.Mirasso
· ¸ ¸ ¸ ¹
a ° ® a ° a ¯
Interpolación y Aproximación
x
0 1 2
½ ° ¾ ° ¿ Año 2014
° ® ° ¯
x0
½ ° ¾ ° ¿
x1
x2
11 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
a1 a2 de donde se tiene que a0 El Polinomio Interpolante usando la Base de Polinomios de Newton es
P2 ( x)
Dr. Ing. A.Mirasso
Interpolación y Aproximación
Año 2014
12 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
FUNCIÓN ERROR DE INTERPOLACIÓN Dados (n+1) puntos de coordenadas (xi; yi=f(xi)) con i=0, n, la Función Error de Interpolación E(x) es la diferencia entre la función f(x) y y el polinomio de interpolación Pn(x), y se puede expresar en la forma: E(x) = f(x) - Pn(x)= f(x) -6 ak Ik(x),
con Ik(x) n funciones bases conocidas (Lagrange, Newton, etc); y ak son los coeficientes que hacen cero el residuo. La función Error de interpolación E(x) tiene (n+1) ceros, en cada uno de los xi datos,
f(x) Pn(x)
por lo que se puede expresar como un polinomio de al menos grado (n+1), como E(x)
E(x)= C (x-x0) (x-x1) (x-x2)…… (x-xn)
La C se determinará de modo que la función auxiliar W(x) sea cero para todo valor de x. W(x)= f(x) - Pn(x)-E(x)
x0
x xi
xn
La constante C se determinará de modo que la igualdad f(x)=Pn(x)+E(x), se cumpla para cualquier valor de x. La función W (x) vale cero en cada xk dato (k=0,n), es decir tiene n+1 ceros. Pero para cualquier otro xi z xk, se elige C de modo que W (xi)=0, entonces d 1W d 2W tiene (n+2-1) ceros, para cada x ; ĺ tiene (n+2-2) ceros, para cada xi; i dx 1 dx 2 d n 1W d n 1 f d n 1 Pn d n 1 E ĺ tiene (n+2-(n+1)) cero, para cada xi. Esta derivada es: n 1 dx n 1 dx n 1 dx n 1 dx
W(x) tiene (n+2) ceros, para cada xi; ĺ
!!. d
n 1
W dx n 1
d n 1W dx n 1
d n 1 f 0 C ( n 1)" dx n 1
Así para cada valor xi z xk existe una C que hace el error de interpolación se pueda expresar E ( x)
d n 1 f ( x ) dx n 1 X
[
1 ( x x0 )( x x1 )( x x2 ).............( x xn ) (n 1)"
con [ H (x0, xn),
Esta expresión establece que: x El error de interpolación en las abscisas datos es cero x La interpolación es exacta si f(x) es un polinomio hasta de grado n. Un tratamiento detallado se puede ver en: Análisis Numérico, R. Burden y D. Faires (1998). Capitulo 3, Teorema 3.3.; Análisis Numérico, W. Smith (1988). Capitulo 7, Teorema 2. Dr. Ing. A.Mirasso
Interpolación y Aproximación
Año 2014
13 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE MÍNIMOS CUADRADOS
Dada una función discreta, y=f(x) de R o R definida mediante (n+1) puntos (xi; yi=f(xi)) con i=0,n. y Se propone P (x)= 6 a I (x) con k=0,m.
ak son las incógnitas; y {Ik(x)} es una Base elegida ri = f(xi) - Pm(xi) con k=0,m; Se define: Residuo m
k
k
ri
Donde
r0 ½ °r ° °° 1 °° ®r2 ¾ °... ° ° ° °¯rn °¿
ri =
yi
y 0 ½ § I 0 ( x0 ) ° y ° ¨ I (x ) °° 1 °° ¨ 0 1 ¨ ® y 2 ¾ ¨ I0 ( x2 ) ° ... ° ¨ .... ° ° ¨ °¯ y n °¿ © I0 ( x n )
r
-
I1 ( x0 ) I1 ( x1 ) I1 ( x 2 ) .... I1 ( x n )
6 ak Ik(xi) I 2 ( x0 ) I 2 ( x1 ) I 2 ( x2 )
.... ....
.... I2 ( xn )
.... ....
y )a
i=0,n
x xi
xn
y
¸ .... ¸° ... ° ° ° I m ( x n ) ¸¹°¯a m °¿
>
La condición para que la Suma de los Cuadrados de los Residuos tome un valor extremo es
I2(xi)
x
@
x0
¦ 2 r (a ) I ( x ) i
j
j
i
xi rT I j
0
xn
j 1, m
i 1,n
Como se debe cumplir para todo j=1,m; finalmente resultan las siguientes ECUACIONES NORMALES
ĭ T ĭa
ĭT y ,
Sistema de ecuaciones lineales cuya solución da los coeficientes ak. Así la APROXIMACIÓN resulta Pm(x)= Dr. Ing. A.Mirasso
Interpolación y Aproximación
Año 2014
I0(xi)
I1(xi)
min 6(r (a ) ) 2 ,
Los coeficientes ak son tales que minimizan la Suma de los Cuadrados de los Residuos 2 i k 2
min r
rn Pm(x)
r2
x0
I m ( x0 ) · a 0 ½ ¸ I m ( x1 ) ¸°° a1 °° ° ° I m ( x 2 ) ¸® a 2 ¾
....
r1
r0
f(x)
6 ak Ik(x)
con k=0,m 14 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE MÍNIMOS CUADRADOS
Dada una función discreta, y=f(x) de R o R definida mediante (2) puntos (xi; yi=f(xi)) con i=0,1. Con la Base
r1 ½ ° ° ®r2 ¾ °r ° ¯ 3¿
= {1, x} el vector residuo r
° ® ° ¯
· ¸ a 0 ½ ¸® ¾ ¸¯ a1 ¿ ¹
½ § ° ¨ ¾¨ ° ¨ ¿ ©
y )a
T ĭ ĭa El sistema de ecuaciones normales
§ ¨¨ ©
·a0 ½ ¸¸® ¾ ¹¯ a1 ¿
EJEMPLO
y=f(x) 3 4
ĭT y resulta
½ ® ¾ ¯ ¿ de donde se tiene que
1
resulta
x 1,5 3
a0
a1
El Polinomio de Aproximación usando la Base de Polinomios elementales {1, es Pa( x) Dr. Ing. A.Mirasso
x}
x Interpolación y Aproximación
Año 2014
15 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTEGRACIÓN NUMÉRICA
Planteo de Problema Integración de Newton Cotes Método de Trapecios Simple y Múltiple Regla de Integración y Error Método de Simpson Simple y Múltiple Regla de Integración y Error Extrapolación de Richardson e Integración de Romberg Integración de Gauss Legendre Ejemplos
Dr. Ing. A.Mirasso
Integración Numérica
Año 2014
1 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTEGRACIÓN NUMÉRICA El propósito es abordar el cálculo de integrales definidas de funciones. Se asume que:
la función y = f(x): R ĺ R, no singular, continua (al menos por tramos), es conocida en forma Discreta Analítica f(x) se conoce en todo x [x0;xn]. f(x) se conoce en xi, con i = 0,1,2,...,n y
y
f(x)
f(x)
x x0
xi
x
xn
x0
xi
xn
Se busca mostrar que:
³ f ( x)dx ¦ w
¦w
es posible calcular la integral definida de la y=f(x) como una combinación lineal: b
I
k
a
k 0, N
f ( xk )
k
k 0,N
yk ,
donde los coeficientes wk son valores particulares para cada regla de integración y los yk=f(xk) son los valores de la función discreta.
Dr. Ing. A.Mirasso
Integración Numérica
Año 2014
2 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTEGRACIÓN NUMÉRICA Si f(x) está dada en forma discreta es posible interpolar f(x) colocando un polinomio Pn(x), de grado n, por los (n+1) puntos datos. Si f(x) está dada en forma analítica se pueden "extraer" esos (n+1) puntos, y tener la versión discreta de f(x). Es posible expresar a f(x) como suma del Polinomio de y Interpolación mas la función Error de Interpolación
f ( n1) ([ ) f(x) = Pn(x) + (n 1)! ( x x0 )( x x1 )........( x xn )
f(x) Pn(x)
Resulta posible obtener la integral en la forma
³ f ( x)dx ³ P ( x) H Xn
Xn
I
n
³
¦w
( x) dx
Xn
f ( x) dx
k
³ P ( x)dx ³ H Xn
n
X0
X0
X0
I
n
Xn
yk E n
( x) dx
E(x)
X0
In En ,
x0
x xi
xn
k 0, N
X0
Resultando, la aproximación de la integral In y su Error de Integración En en la forma:
³ P ( x)dx ¦ w
Xn
In
n
n
X0
Dr. Ing. A.Mirasso
k
yk
k 0, N
³H
En
n
( x) dx
X0
X0
Integración Numérica
³
Xn
Xn
Año 2014
f ( n 1) ([ ) ( x x0 ) ( x x n ) dx (n 1)!
3 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTEGRACIÓN NUMÉRICA DE NEWTON COTES Las reglas de integración de Newton Cotes se basan en interpolar con Polinomios de LAGRANGE. Para los n+1 puntos datos, el polinomio interpolante es y n
¦ y l ( x)
Pn ( x)
i
f(x)
i
Pn(x)
i 0
³ ¦ y l ( x)dx ¦ ³ y
Así el valor aproximado de la integral Xn
Xn n
In
i
¦
X0 i 0
In
i
yk ³ lk ( x) dx
¦y
lk ( x) dx
k 0, N X 0
Xn
k 0, N
k
X0
resulta
³l
Xn
wk
k
wk
k 0, N
( x) dx
X0
k
³ (x
Xn n
X0 j 0 jzk
Así el Error de la aproximación de la integral
En
³H
n
( x) dx
X0
Dr. Ing. A.Mirasso
(x x j ) k
Integración Numérica
lk(x)
xj)
f ( n 1) ([ ) n ( x x0 ) ( x x n ) dx (n 1)! X³0 X
Xn
l0(x)
dx
x0
x0
I0(xi) xi
x xn
Ik(xi) x
x0
xn
xi
f ( n 1) ([ ) n 2 h D n 1 (n 1)! Año 2014
4 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
REGLA DE INTEGRACIÓN DE LOS TRAPECIOS
Es una regla de integración de Newton ! Cotes por 2 puntos (xi;yi), (xi+1;yi+1). La integral
³ f ( x) dx
xi 1
I
f(x)
xi
³ >P ( x) H
( x ) @dx
³ >y
l i ( x ) y i 1 l i 1 ( x ) @dx
Se resuelve con un polinomio interpolante degrado uno xi 1
I
1
xi
donde
li ( x)
x x i 1 x i x i 1
li 1 ( x)
x xi xi 1 xi
³
xi 1
I
I
yi
³
xi 1
xi
1
xi 1
x x i 1 x i 1 x i
xi
1
( x ) dx ,
xi
Xi
Xi+1
l0(x)
li(x)
li+1(x)
es una recta que vale 0 en xi y 1 en xi+1.
ª x xi º x xi 1 y y i 1 » dx « i xi 1 xi ¼ xi xi 1 ¬
³
xi 1
xi
x xi dx xi 1 xi
I
³ H ( x)dx,
xi 1
³ H ( x)dx xi
Xi
Xi+1
yi wi yi 1 wi 1 E1
1
xi
y º ªy hi « i i 1 » E1 2 ¼ ¬2
Integración Numérica
X
1
xi 1
y, llamando paso hi a la diferencia xi+1-xi y operando, se puede llegar a
Dr. Ing. A.Mirasso
³H
X
es una recta que vale 1 en xi y 0 en xi+1,
x xi 1 dx yi 1 xi xi 1 xi
i
xi 1
I 1 E1 ,
Año 2014
5 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
ERROR DE REGLA DE INTEGRACIÓN DE LOS TRAPECIOS
E1
xi 1
1
[ ( xi ; xi 1 ). xa
siendo
f ( 2 ) ([ ) ( x xi )( x xi 1 ) dx 2!
³ H ( x)dx ³
El error de interpolación es xb
xi
x xi xi 1 xi
t ( 1) 1 ( 1)
f ( 2 ) ([ ) 2!
³ ( x x )( x x
xi 1
i
i 1
) dx f(x)
xi
P1(x)
Se propone hacer un cambio de variable mediante
x xi
(x-xi) (x-xi+1) X
,
Xi
h (t 1) / 2
Definiendo h=xb-xa, es posible despejar,
x xi 1 dx
Xi+1
h(t 1) / 2
t
( h / 2) dt
Así se puede hallar
³ ( x xi )( x xi1 )dx
xi 1
xi
h(t 1) h(t 1) h ³1 2 2 2 dt 1
h 8
2 ³ (t 1)dt
3 1
1
-1 h 3
1
1 6
de manera que, sustituyéndola en (5), se obtiene
E1 Dr. Ing. A.Mirasso
f ( 2 ) ([) § 1 · 3 ¨ ¸h 2! © 6 ¹ Integración Numérica
1 3 ( 2) h f ([) 12
para cierto punto [ (xa, xb).
Año 2014
6 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
REGLA DE INTEGRACIÓN DE LOS TRAPECIOS MÚLTIPLES
Se busca I R,
³ f ( x)dx .
y
Xn
I
f(x)
X0
P1(x)
Se divide el intervalo [x0; xn] en subintervalos [x0; x1], [x1; x2], }, [xn-1; xn]. Así
³ f ( x) dx ³ f ( x)dx ³ f ( x) dx .
X1
I
Xn
X2
x x0
X n 1
xn En cada uno de los n subintervalos se aplica la regla de los trapecios, se aplica trapecios simple X0
I
h0
y 0 y1 2
X1
E1 (h0 ) h1 3
y1 y 2 2
E1 (h1 ) .............. hn 1 3
y n1 y n 2
E1 (hn 1 ) . 3
Si todos los intervalos tienen igual longitud hi=h, esa fórmula se simplifica y se tiene la regla de trapecios múltiple:
I
y º ª y 0 n 1 h « ¦ y i n » E 1M 2¼ ¬2 i1
n 1 hª º y 2 y y ¦ 0 i n » E 1M , « 2¬ i 1 ¼
donde E1M es el error total que se acumula al sumar los n errores provenientes de la aplicación de la regla en cada subintervalo, y está dado por
E1M Dr. Ing. A.Mirasso
Integración Numérica
xn x0 h 2 f ( 2) ([ ) 12
Año 2014
7 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
EJEMPLO DE LA REGLA DE INTEGRACIÓN DE LOS TRAPECIOS MÚLTIPLES
³ sen(x)dx cuyo resultado exacto es I = -cos x S
Se busca resolver
I
n 1 h ª º I ³ sen( x )dx y y y 2 ¦ n» i §2 « 0 i 1 ¼ ¬ 0
0
S
donde X 0 ʌ/4 ʌ/2 3ʌ/4
S 0
= cos 0 cos S = 2
h S , 2
yi = sen (xi), i=0,},n. Y=sen(x) h1 = ʌ h2 = ʌ/2 h3 = ʌ/4 0 1 1 1 0 0 2 2/2
ʌ
1 2/2
0 0
2 0
2 2
0 0
1 S1
1 S2
1 S3
Es fácil concluir que la mejor aproximación es I3. Los errores cometidos por cada aplicación de la regla son, respectivamente O(h13), O(h22) y O(h32). Dr. Ing. A.Mirasso
Integración Numérica
Año 2014
8 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
³0 sen(x)dx S
Extrapolación de Richardson para el ejemplo
I
A partir de dos valores aproximados de la Integral, obtenidos con la misma Regla y pasos distintos, es posible encontrar un valor mejorado mediante Extrapolación de Richardson.
I Dr. Ing. A.Mirasso
E I ( h 2 ) I( h 1 ) E 1
Integración Numérica
con
Año 2014
E
§ h1 ¨¨ © h2
· ¸¸ ¹
n
9 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
³ sin(2x) dx
3S / 2
Ejemplo 2 Regla de Integración de los Trapecios Múltiples Int
0
n h w 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 Integral h
16 0,294524311 w* Sin(2x) 0 1,111140466 1,847759065 1,961570561 1,414213562 0,390180644 -0,765366865 -1,662939225 -2 -1,662939225 -0,765366865 0,390180644 1,414213562 1,961570561 1,847759065 1,111140466 2,1439E-15 0,970916536 0,294524311
n 8 h 0,58905 w w* Sin(2x) 1 0
n 4 h 1,1781 w w* Sin(2x) 1 0
2
1,84776
2
1,41421 2
2
-0,7654
2
-0,7654
2
1,41421 2
2
1,84776
1
2,1E-15 1 0,88157 0,58905
3S / 2
n 1 h 4,7124 w w* Sin(2x) 1 0
-2 2
-2
1,4142
2E-15 1 0,488 1,1781
2,1E-15 1 -2,3562 2,35619
2E-15 5E-15 4,7124
1,50000
1,50000
1,50000
1,50000
1,00000
1,00000
1,00000
1,00000
0,50000
0,50000
0,50000
0,50000
0,00000
0,00000
0,00000 0
1
2
3
4
-0,50000 seno(2*x) -1,00000 -1,50000
Dr. Ing. A.Mirasso
5
0
1
2
3
4
-1,00000
Forma discreta
0,00000 0
5
1
seno(2*x)
-1,00000
Integración Numérica
3
4
-1,50000
Año 2014
5
0
1
2
3
4
-0,50000 seno(2*x) Forma discreta
Forma discreta -1,50000
2
-0,50000
-0,50000
1
1,4142
-2 2
2
n 2 h 2,35619 w w* Sin(2x) 1 0
1 cos(2 x) 2 0
-1,00000
seno(2*x) Forma discreta
-1,50000
10 de 15
5
Facultad de Ingeniería, Universidad Nacional de Cuyo
³ sin(2x) dx
3S / 2
Extrapolación de Richardson para el ejemplo Int
E I( h 2 ) I( h 1 ) E 1
I
0
con
1 cos(2 x) 2 0
3S / 2
§ h1 ¨¨ © h2
E
· ¸¸ ¹
1
n
Y la extrapolación de Richardson sucesiva (Método de Romberg) resulta: n
Orden
h^2
h^4
h^6
h^8
Intervalos
Paso h
I(h^2)
I(h^4)
I(h^6)
I(h^8)
I(h^10)
16
0,29452
0,97092
8
0,58905
0,88157
1,00069753
4
1,1781
0,48798
1,01277013
0,999892687
2
2,35619
-2,3562
1,43604331
0,98455192
1,000136191
1
4,71239
5,1E-15
-3,14159265
1,741219036
0,972541331
1,50000
1,50000
1,50000
1,50000
1,00000
1,00000
1,00000
1,00000
0,50000
0,50000
0,50000
0,50000
0,00000
0,00000
0,00000 0
1
2
3
4
-0,50000 seno(2*x) -1,00000 -1,50000
Dr. Ing. A.Mirasso
5
0
1
2
3
4
-1,00000
Forma discreta
0,00000 0
5
1
seno(2*x)
-1,00000
Integración Numérica
3
4
-1,50000
Año 2014
5
0
1
2
3
4
-0,50000 seno(2*x) Forma discreta
Forma discreta -1,50000
2
-0,50000
-0,50000
1,00024441
-1,00000
seno(2*x) Forma discreta
-1,50000
11 de 15
5
Facultad de Ingeniería, Universidad Nacional de Cuyo
REGLA DE INTEGRACIÓN DE SIMPSON Es una cuadratura de Newton ! Cotes con n = 2, es decir con tres puntos. Se interpola mediante un polinomio de Lagrange de grado dos y luego se integra en forma aproximada ese polinomio.
³ f ( x)dx ,
X2
I
f(x)
X0
Si
f(x) =¦Yi li(x) + E2(x),
H 2 ( x)
f
x x 0 x x1 x x 2
( 3)
¦yZ 3"
2
I2
x x0 x x2 l 2 ( x) x1 x 0 x1 x 2
x x1 x x 2 l1 ( x ) x 0 x1 x 0 x 2
l 0 ( x)
i
i 0
con i = 0,1,2,
i
Zi
³ l ( x)dx, i
x2
i
³H
x x 0 x x1 x 2 x 0 x 2 x1
x0
X0
X1
X2
l0(x)
X X0
x2
0,1,2. E 2
X
2
( x) dx
X1
X2
l1(x)
x0
Entonces, si los intervalos son iguales (h1= h2 = h), se tiene:
I
4 1 º h5 ( 4) ª1 h « y0 y1 y2 » f ([ ) 3 3 ¼ 90 ¬3
Dr. Ing. A.Mirasso
Integración Numérica
X X0
X1
X2
l2(x)
X X0
Año 2014
X1
X2
12 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
INTEGRACIÓN DE GAUSS LEGENDRE
³
h ³ G (t ) dt
Se busca resolver la b
I
1
f ( x) dx
a
G (t ) f (c (
1
ba ) t) 2
n § n · h¨ ¦ Z i G (t i ) R ¸ h¦ Z i G (t i ) E i 0 ©i 0 ¹ ba x c( )t 2
El problema se resuelve en el dominio unitario >-1, 1@, mediante Mapeo o Cambio de variables conveniente.
³ G(t)dt 1
1
Z1G ( t 1 ) Z 2 G ( t 2 ) R
Z1, Z2, t1, t2 deben ser tales que el error de integración sea R = 0 para polinomios de hasta 3º grado. Comparando los resultados exactos de la integral con el resultado propuesto por la regla, se tiene:
³ 1dt 1
1
³ tdt
2 1Z1 1Z 2
1
³ 1
1
t 2 dt
³ t dt
0 2 3
1
3
1
Z1
t1Z1 t 2 Z 2
1
0
t1 Z1 t 2 Z 2 2
2
t1
Z2 t2
1
I 3 3
a
c
b
ª § 3· § 3 ·º ¸» E ¨ ¸ h «G ¨ G¨¨ ¸ ¸ «¬ © 3 ¹ © 3 ¹»¼
t13Z1 t 2 3Z 2
Dr. Ing. A.Mirasso
Integración Numérica
Año 2014
13 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
Integración de Gauss Legendre. Generalización. Se busca resolver la
³ G(t )dt ¦ Z G(t ) R 1
n
1
i 0
i
i
En una regla de n+1 puntos existen 2n+2 coeficientes a determinar, que son los (n+1) wk y las (n+1) abscisas tk. Se exige que la regla de integración sea exacta cuando integra polinomios de hasta grado 2n+1. Dada una función polinómica G(t) de hasta grado (2n+1), existen polinomios Cn(t) y rn(t) de grado menor o igual a n, tales que G(t)= Cn(t) pn+1(t )+ rn(t) Siendo pn+1(t ) un polinomio de grado (n+1), que en particular puede ser un polinomio de Legendre. Así es posible expresar la integral
³ G (t ) dt ³ (C 1
1
I
1
1
n
³C
(t ) p n 1 (t ) rn (t )) dt
1
1
(t ) p n 1 (t )dt ³ rn (t ) dt 1
n
1
³ r (t ) dt ¦ w 1
1
n
k
rn (t k )
k 0, n
Ya que cualquier polinomio de grado n es ortogonal al polinomio de Legendre de grado n+1 en el dominio [-1,1]. Para elegir las abscisas tk, es posible elegir las (n+1) raíces del polinomio de Legendre pn+1(t) de grado (n+1), que se las denomina tr, y en las que se verifica que G(tr)= Cn(tr) pn+1(tr )+ rn(tr)= rn(tr) Así la integral resulta
³ G (t ) dt ³ r (t ) dt ¦ w 1
I
1
Dr. Ing. A.Mirasso
1
1
n
k
k 0, n
rn (t k )
¦w
k
G (t k )
k 0, n
Integración Numérica
Año 2014
14 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
Para calcular los wk, Es posible plantear una interpolación del polinomio de grado n rn(t), con polinomios de Lagrange, tomando como abscisas conocidos las (n+1) raíces tk del polinomio de Legendre de grado (n+1) rn(t)= rn(t0) l0(t) + rn(t1) l1(t) + rn(t2) l2(t) +.... rn(tn) ln(t) o bien considerando que G(tr)= rn(tr) rn(t)= G(t0) l0(t) + G(t1) l1(t) + G(t2) l2(t) +.... G(tn) ln(t)
I
³ G (t ) dt ³ r (t ) dt 1
1
1
1
n
1 ª º G ( t ) dt G ( t ) l ( t ) dt ¦« k ³ k » ³ k 0, n ¬ 1 1 ¼ 1
I
º ª G ( t ) l ( t ) « ³1 ¬k¦0, n k k »¼dt 1
¦ wk G(tk )
º ª1 G ( t ) l ( t ) dt » « ¦ k k ³ k 0 , n ¬ 1 ¼
³ l (t ) dt 1
Con
wk
1
k 0, n
k
Es decir que los coeficientes wk son las integrales de los polinomios de Lagrange de grado n+1, definidos por las (n+1) raíces tk de los polinomios de Legendre de grado n+1. Así resulta,
³ G (t ) dt ¦ w
Siendo
1
I
1
k
G (t k )
k 0, n
³ l (t ) dt
wk
1
k
Se toma (n+1) puntos de Gauss (raices del polinomio de grado n+1), y se puede integrar en forma exacta hasta polinomios de grado (2n+1) (n+1)Puntos de Gauss 1 2 Grado exacto 1 3
Dr. Ing. A.Mirasso
Integración Numérica
Y (tr )raíces del Polinomio de Legendre pn+1(tr )=0
1
Año 2014
3 5
4 7
15 de 15
Facultad de Ingeniería, Universidad Nacional de Cuyo
DERIVACIÓN NUMÉRICA
Planteo de Problema Derivadas a partir de Interpolación con Polinomios de Newton Derivada Primera Derivada Segunda Derivadas a partir de Serie de Taylor Derivada Primera Derivada Segunda Extrapolación de Richardson Ejemplos
Dr. Ing. A.Mirasso
Derivación Numérica
Año 2014
1 de 8
Facultad de Ingeniería, Universidad Nacional de Cuyo
DERIVACIÓN NUMÉRICA El propósito es abordar el cálculo de derivadas de distinto orden de funciones discretas. Supuestos
La función y = f(x): R ĺ R, no singular, continua, es conocida en forma Discreta Analítica f(x) se conoce en xi, con i = 0,1,2,...,n f(x) se conoce en todo x [x0;xn]. y f(x)
y f(x)
x x0
xi
x
xn
xi
x0
xn
Hipótesis
¦ ck f ( xk )
¦ ck y k
Es posible calcular la derivada n-ésima de la y=f(x) evaluada en Xj como una suma:
D
d n ( f ( x )) dx n
Xj
k 0, N
k 0, N
& & cT y
,
Donde los coeficientes ck son valores particulares para cada regla de derivación y los yk=f(xk) son los valores de la función discreta. Dr. Ing. A.Mirasso
Derivación Numérica
Año 2014
2 de 8
Facultad de Ingeniería, Universidad Nacional de Cuyo
DERIVACIÓN NUMÉRICA Si f(x) está dada en forma discreta es posible interpolar f(x) colocando un polinomio Pn(x), de grado n, por los (n+1) puntos datos. Si f(x) está dada en forma analítica se pueden "extraer" esos (n+1) puntos, para la versión discreta de la función f(x).
f ( n1) ([ ) f(x) = Pn(x) + (n 1)! ( x x0 )( x x1 )........( x xn )
Es posible expresar:
suma del Polinomio de Interpolación mas la función Error de Interpolación. y f(x) Pn(x) E(x)
D
d n ( f ( x)) dx n
Dr. Ing. A.Mirasso
¦ ck f ( xk ) E n x0
Resulta posible
Xj
k 0, N
Derivación Numérica
x xi
xn
¦ ck y k E n
& & c T y En
k 0, N
Año 2014
3 de 8
Facultad de Ingeniería, Universidad Nacional de Cuyo
DERIVACIÓN NUMÉRICA
Si f(x) está dada en forma discreta con (n+1) puntos, es posible interpolar f(x) colocando un polinomio de grado n. A partir de
f(x) = Pn(x) +
f ( n1) ([ ) ( x x0 )( x x1 )........( x xn ) (n 1)!
d n Pn ( x ) H n ( x ) dx n Xj
Es posible obtener derivadas hasta de grado n, y evaluarlas en cualquier punto Xj, de la forma: D
D
d n ( f ( x )) dx n d
n
Pn ( x ) dx n
D
Xj
Xj
Dn
Pn ( x )
Resultando Dn
d
n
dx
n
Xj
y
d H n ( x ) dx n Xj
f(x)
n
Pn(x)
En
¦c
k 0,N
k
f ( xk ) k
¦c
k
yk
& & cT y
0,N
con el Error de Derivación dado por E n
d n H n ( x ) dx n
E(x)
x0
x xi
· § f ( n 1) ([ ) d ¨¨ ( x x 0 ) ( x x n ) ¸¸ ¹ © ( n 1)! n dx
xn
n
Xj
Pregunta: ¿Qué cantidad de puntos se necesita para evaluar una derivada primera?, ¿y para una derivada segunda? Dr. Ing. A.Mirasso
Derivación Numérica
Año 2014
4 de 8
Facultad de Ingeniería, Universidad Nacional de Cuyo
DERIVADA PRIMERA
Si f(x) está dada en forma discreta con 2 puntos (x0;y0), (x1;y1) es posible interpolar f(x) colocando un polinomio de grado 1, usando polinomios de Newton, en la forma: f ( 2 ) ([ ) ( x x0 ) ( x x1 ) f ( x ) a0 a1 ( x x0 ) y ( 2)! La derivada primera de f(x) es d ( f ( x )) f ( 2 ) ([ ) d (( x x0 ) ( x x1 )) a1 dx ( 2)! dx d ( f ( x )) dx
f ( 2 ) ([ ) a1 [( x x1 ) ( x x0 )] ( 2)!
P1(x) f(x)
x y
x0
x1
n0(x)
1
Cuando se evalúa la derivada primera de f(x) en x0 y en x1 se obtienen: Derivada Primera Adelante y1 y 0 f ( 2 ) ([ ) d ( f ( x )) 'x con a1 a1 y 'x x1 x0 ( 2)! 'x dx X X0 Derivada Primera Atrás y1 y 0 f ( 2 ) ([ ) d ( f ( x )) 'x con a1 y 'x x1 x0 a1 'x ( 2)! dx X X1
n1(x) x
x1
x0
Derivada Primera es constante en el intervalo; y el Error es lineal. Dr. Ing. A.Mirasso
Derivación Numérica
Año 2014
5 de 8
Facultad de Ingeniería, Universidad Nacional de Cuyo
DERIVADA SEGUNDA
Si f(x) está dada en forma discreta con 3 puntos (x0;y0), (x1;y1) es posible interpolar f(x) colocando un polinomio de grado 2, usando polinomios de Newton, en la forma: f ( 3) ([ ) f ( x ) a0 a1 ( x x0 ) a 2 ( x x0 ) ( x x1 ) ( x x0 ) ( x x1 ) ( x x2 ) (3)! y La derivada segunda de f(x) es f(x) d 2 ( f ( x )) f ( 3) ([ ) d 2 (( x x0 ) ( x x1 ) ( x x2 )) 2a 2 2 Pn(x) (3)! dx dx 2 y 2 y1 y1 y 0 x x 2 x1 x1 x 0 1 y 2 y1 y1 y 0 xn 2a 2 2 2 Con n (x) 2 'x 'x x2 x0 2'x n1(x) y Con lo que la derivadas segunda es n0(x) d 2 ( f ( x )) 1 ( y 0 2 y1 y 2 ) E D ( x ) 2 2 dx 'x x
E D ( x)
f ( 3) ([ ) d 2 (( x x0 ) ( x x1 ) ( x x2 )) (3)! dx 2
x0
x1
x2
Derivada Primera es constante en el intervalo; y el Error es lineal. ¿El error será igual a cero en alguno de los puntos datos? Dr. Ing. A.Mirasso
Derivación Numérica
Año 2014
6 de 8
Facultad de Ingeniería, Universidad Nacional de Cuyo
DERIVADA PRIMERA EN BASE A SERIE DE TAYLOR
Si f(x) está dada en forma discreta con 2 puntos (xs;ys), (xs+1;y s+1) es posible considerar el desarrollo de Serie de Taylor de y s+1=f(x s+1) alrededor de xs P1(x) en la forma:
f s 1
h2 h3 h 4 cv f s h f sc f scc f sccc f s O h5 2! 3! 4!
y
f(x) ,
de donde se puede despejar la derivada primera en xs
ys+1
ys xs
2 3 4 nh f cc xs nh f ccc xs nh f xs nh f xs n h f c xs
2!
Dr. Ing. A.Mirasso
Derivación Numérica
3!
Año 2014
x
xs+1
f cv xs O h5 4!
7 de 8
Facultad de Ingeniería, Universidad Nacional de Cuyo
DERIVADA SEGUNDA EN BASE A SERIE DE TAYLOR
Si f(x) está dada en forma discreta con 3 puntos (xs-1;y s-1) (xs;ys), (xs+1;y s+1) es posible considerar el desarrollo de Serie de Taylor de y s±1=f(x s±1) alrededor de xs en la forma: y
f s 1 f s 1
h2 h3 h 4 cv f s h f sc f scc f sccc f s ...... , 2! 3! 4!
f(x)
h2 h3 h 4 cv f s h f sc f scc f sccc f s ..... 2! 3! 4!
P2(x)
>D f s 1 E f s J f s 1 @
x
y con ello se puede plantear en xs la derivada segunda de f(x) como
fs ''
2 3 4 nh f cc xs nh f ccc xs nh f xs nh f xs n h f c xs
2!
Dr. Ing. A.Mirasso
Derivación Numérica
3!
Año 2014
xs-1
xs
xs+1
f cv xs O h5 4!
8 de 8
Facultad de Ingeniería, Universidad Nacional de Cuyo
SOLUCIÓN NUMÉRCIA CONTORNO
DE
ECUACIONES
DIFERENCIALES
CON
VALORES
Ecuaciones diferenciales: una clasificación Ecuaciones Diferenciales Ordinarias Orden Superior Problema de Valores de Contorno Obtención de Sistemas de Ecuaciones Lineales No Homogeneos Obtención de Sistemas de Valores y vectores propios Problema de Valores Iniciales Ecuaciones Diferenciales a Derivadas Parciales
Ejemplos Resumen
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
1 de 14
DE
Facultad de Ingeniería, Universidad Nacional de Cuyo
SOLUCIÓN NUMÉRCIA DE ECUACIONES DIFERENCIALES ORDINARIAS Se busca obtener una solución aproximada en forma discreta Clasificación de EDO
Primer orden Siempre son de valor inicial
du (t ) A u (t ) 0, dt u (t0 ) U 0 ,
u(t)
A R,
u0
u(t) U1
La solución discreta es tal que aproxima a la solución exacta del problema
Uk # u(tk)
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
t0
t1
Año 2014
U2 t2
U3
t
t3
2 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
SOLUCIÓN NUMÉRCIA DE ECUACIONES DIFERENCIALES ORDINARIAS Se busca obtener una solución aproximada en forma discreta Orden Superior Pueden ser de valores iniciales de valores de contorno
d 2 u (t ) g u (t ) 0, 2 L dt du (t ) u (t 0 ) u 0 , dt t 0
d 2 u ( x) u ( x) x 2 dx u (0) 0
v0
u(t)
0
en :
^x H R : 0 d x d 1` u (1)
0
U(x)
u0 U1
U2
t1
t2
U3
t
U0
t3
Uk # u(tk)
U1
U2
X1
X2
U3 X=3
U4
X
X=1
Uk # u(xk)
La solución discreta es tal que aproxima a la solución exacta del problema
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
3 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
DISCRETIZACIÓN DE ECUACIONES DIFERENCIALES ORDINARIAS
Se busca u(x) solución de d 2 u ( x) u ( x) x 2 dx u ( 0) 0 u (1)
0
en :
^x H R : 0 d x d 1`
0
se plantea encontrar una solución aproximada en forma de función discreta, sólo en algunos puntos elegidos del dominio : , equidistantes entre sí, identificados con su abscisa Xk. Es decir, se busca U(Xk)=Uk con k=0,N; función discreta que es una aproximación de la función continua u(x).
U(x)
U0
U1 X=0,25
Dr. Ing. A.Mirasso
U2 X=0,5
U3 X=0,75
Ec.Dif. con Valores de Contorno
U4
X
X=1
Año 2014
4 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
Se busca la función discreta U(Xk)=Uk con k=0,N; que es una aproximación de la función continua u(x). d 2 u ( x) u ( x) x dx 2 u (0) 0 u (1) 0
0
en :
^x H R : 0 d x d 1`
U(x)
U1
U0
X=0,25
U2 X=0,5
U3 X=0,75
U4
X
X=1
Se divide el dominio : en N segmentos iguales En cada uno de los Xk se puede plantear la ecuación diferencial a resolver pero con una aproximación de la derivada segunda en forma de derivada numérica
1 >U k 1 2 U k U k 1 @ U k X k 'x 2
0
en X k
con k
1, ( N 1)
De estas ecuaciones se pueden plantear tantas como puntos interiores; es decir N-1 ecuaciones y se tienen N+1 incógnitas. Además se tienen las dos ecuaciones correspondientes a las Condiciones de Contorno, que agregan dos ecuaciones más. Así se tienen N+1 ecuaciones con N+1 incógnitas.
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
5 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
d 2 u ( x) u ( x) x dx 2 u (0) 0 u (1) 0
0
en :
^x H R : 0 d x d 1`
U(x)
U1
U0
X=0,25
U2 X=0,5
U3 X=0,75
U4
X
X=1
Si se consideran 5 puntos en todo el dominio, que son 3 en el interior y uno en cada uno de los bordes, es posible plantear U0
en X0 se debe cumplir que: en X1 se debe cumplir que:
en X2 se debe cumplir que:
en X3 se debe cumplir que:
1 >U 0 2 U1 U 2 @ U1 X 1 0,252
0
1 >U 2 2 U 3 U 4 @ U 3 X 3 0,25 2 U4
0 · §1 0 ¸ ¨ 1¸ ¨ 0 1 2 ¸¹ ¨© 0 0 0,03488525½ ° ° ®0,05632582¾ °0,05003676° ¿ ¯
O bien en forma de sistema de ecuaciones lineales
Dr. Ing. A.Mirasso
0
1 >U1 2 U 2 U 3 @ U 2 X 2 0,252
en X4 se debe cumplir que:
ª § 2 1 « 1 ¨ « 0,252 ¨ 1 2 ¨ 0 1 «¬ © U 1 ½ ° ° La solución aproximada resulta ®U 2 ¾ °U ° ¯ 3¿
0
0 ·º U1 ½ ¸» ° ° 0 ¸» ®U 2 ¾ 1 ¸¹»¼ °¯U 3 °¿
Ec.Dif. con Valores de Contorno
0
0
0 º U1 ½ ª32 1 16 « 16 32 1 16 » °U ° » ® 2¾ « «¬ 0 16 32 1»¼ °¯U 3 °¿
Año 2014
0,25½ ° ° ® 0,5 ¾ °0,75° ¿ ¯
6 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
d 2 u ( x) u ( x) x dx 2 u (0) 0 u (1) 0
0
en :
^x H R : 0 d x d 1`
U(x)
U0
U1 X=0,25
U 0 ½ °U ° °° 1 °° ®U 2 ¾ °U ° ° 3° °¯U 4 °¿
U2 X=0,5
U3 X=0,75
U4
X
X=1
0 ½ °0,03488525 ° °° °° ®0,05632582 ¾ °0,05003676 ° ° ° °¿ °¯ 0
La solución aproximada completa, que incluye las condiciones de borde es:
U 0(1) ½ ° (1) ° °°U 1 °° (1) ®U 2 ¾ °U (1) ° ° 3 ° °¯U 4(1) °¿
ª 3 4 1 0 « 1 0 1 0 « 1 « 0 1 0 1 (2 * 0,25) « 0 1 0 «0 «¬ 0 0 1 4
0º 0 ½ ° » 0» °0,03488525 °° ° ° 0» ®0,05632582 ¾ » 1» °0,05003676 ° ° ° °¿ 3»¼ °¯ 0
Justificar que con la solución obtenida una aproximación de la función derivada primera de esta solución se puede calcular mediante:
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
7 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
Caso N=2 1 >U 0 2 U 1 U 2 @ U 1 X 1 0,5 2
U0
U1
0
0
en X 0
0
U(x)
en X 1
en X 1
U1
U0
O bien
9 U 1 0,5
U1
0X=0,5
X
X=1
La solución aproximada es U1
1 / 18
0,055555
Así el error respecto de la solución exacta en ese punto es E2 E (abs ) 2
Dr. Ing. A.Mirasso
1 senh(0,5) ) (0,5 18 senh(1)
0,001035002
senh(0,5) senh(0,5) 1 ) /(0,5 (0,5 ) 1,83% senh(1) senh(1) 18
Ec.Dif. con Valores de Contorno
Año 2014
8 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
Caso N=4 1 >U 0 2 U 1 U 2 @ U 1 X 1 0 0,5 2 1 >U 1 2 U 2 U 3 @ U 2 X 2 0 0,5 2 1 >U 2 2 U 3 U 4 @ U 3 X 3 0 0,5 2
U0
U4
0
0
en X 0
en X 1
0,25
en X 2
0,5
en X 3
U(x)
U0 0,75
U1 X=0,2
U2
U3
X=0,5
X=0,7
U4
X
X=1
en X 4
O bien
0 º U 1 ½ ª 33 16 « 16 33 16» °U ° « » ® 2¾ «¬ 0 16 33 »¼ °¯U 3 °¿
La solución aproximada es
U 1 ½ ° ° ®U 2 ¾ °U ° ¯ 3¿
0,25½ ° ° ® 0,5 ¾ °0,75° ¯ ¿
0,03488525½ ° ° ®0,05632582¾ °0,05003676° ¯ ¿
Así el error respecto de la solución exacta en ese punto es E4 E (abs ) 4
Dr. Ing. A.Mirasso
U 2 (0,5
U 2 (0,5
Ec.Dif. con Valores de Contorno
senh(0,5) ) senh(1)
0,000264735
senh(0,5) senh(0,5) ) ) /(0,5 senh(1) senh(1)
Año 2014
0,47%
9 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
Caso N=8 U0
0 en X 0 0 1 >U 0 2 U 1 U 2 @ U 1 X 1 (1 / 8) 2
0
en X 1 1 / 8
1 >U 1 2 U 2 U 3 @ U 2 X 2 0 en X 2 2 / 8 (1 / 8) 2 1 >U 2 2 U 3 U 4 @ U 3 X 3 0 en X 3 3 / 8 (1 / 8) 2
………. U8
0
en X 8
1
O bien
0 0 0 0 0 ª 129 64 « 64 129 64 0 0 0 0 « « 0 64 129 64 0 0 0 « 64 129 64 0 0 0 « 0 « 0 64 129 64 0 0 0 « 64 129 64 0 0 0 « 0 « 0 64 129 0 0 0 0 « ¬«
º U1 ½ » °U ° » ° 2° » °U 3 ° » ° ° » °U 4 ° » ®U 5 ¾ » ° ° » °U 6 ° » °U ° » ° 7° »¼ °¯ °¿
1 / 8 ½ °2 / 8° ° ° °3 / 8 ° ° ° °4 / 8° ® ¾ °5 / 8 ° °6 / 8 ° ° ° °7 / 8° ° ° ¯ ¿
^0,0183367; 0,03500678; 0,04831759; 0,05652399; 0,05780107; 0,05021568; 0,03169615`
La solución aproximada es UT
Así el error respecto de la solución exacta en ese punto es E8 E (abs ) 8
Dr. Ing. A.Mirasso
U 4 (0,5
U 4 (0,5
Ec.Dif. con Valores de Contorno
senh(0,5) ) senh(1)
6,65711E - 05
senh(0,5) senh(0,5) ) ) /(0,5 senh(1) senh(1) Año 2014
0,12%
10 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
Evaluación del Error Al considerar el error para distintos niveles de disicretización; es decir, distinto número de segmentos en que se divide el dominio, se tiene senh(0,5) senh(0,5) senh(0,5) (0,5 ) y E (abs ) (0,5 ) ) /(0,5 E U U N
N /2
N
senh(1)
N /2
senh(1)
senh(1)
Cuyas evaluaciones se presentan en la siguiente Tabla N 'x Uaprox(0,5) EN E(abs)N 2 0,5 0,05555556 0,001035002 1,83% 4 0,25 0,05632582 0,000264735 0,47% 8 0,125 0,05652399 6,65711E-05 0,12% 16 0,0625 0,05657389 1,66672E-05 0,03% Si se asume una relación exponencial entre E(abs) N y 'x, la aproximación por mínimos cuadrados da:
E (abs ) N
C 'x P
e 2, 6168 'x 1,9861
Que indica una relación del orden de 'x # 'x que es el error de truncamiento local de la aproximación de derivada segunda utilizado. 1, 9861
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
2
Año 2014
11 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
DISCRETIZACIÓN DE ECUACIONES DIFERENCIALES A DERIVADAS PARCIALES
Se busca u(x,t) solución de
2 w 2 u ( x, t ) 2 w u ( x, t ) x 12 wt 2 wx 2 u (0, t ) 0
u (1, t ) 0
0 en :
^x H R :0 d x d 1`
u ( x,0) sin(S x) wu ( x,0) 3 wt
Se plantea encontrar una solución aproximada en forma de función discreta en la variable x, aunque continua en la variable t. Se pretende encontrar las funciones Uk(t)=u(Xk,t) con k=0,N, en N+1 puntos elegidos del dominio x, identificados con su abscisa Xk. u(x,t)
Cuando se consideran derivadas parciales respecto a la variable x, para t constante, la función a derivar es una función discreta y se puede hacer derivadas numéricas. Cuando se consideran derivadas parciales respecto a la variable t, para x constante, la función a derivar es una función continua y se puede hacer derivadas analíticas.
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
U1(t) U0(t)
0,25
U2(t)
0,5
U3(t) 0,75
U4(t)
X
X=1
t
Año 2014
12 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
Se busca u(x,t) solución de 2 w 2 u ( x, t ) 2 w u ( x, t ) 12 x wx 2 wt 2 u (0, t ) 0
0 en :
^x H R :0 d x d 1`
u ( x,0) sin(S x) wu ( x,0) 3 wt
u (1, t ) 0
u(x,t) U1(t) U0(t)
0,25
U2(t)
0,5
U3(t) 0,75
U4(t)
X
X=1
En cada uno de los Xk se puede plantear la ecuación diferencial a resolver pero con una aproximación de la derivada segunda t respecto a x en forma de derivada numérica; y con la derivada respecto de la variable t en forma analítica evaluada en esa abscisa Xk. Así se puede escribir: U 0 (t ) 0 en X0 se debe cumplir que: en X1 se debe cumplir que: en X2 se debe cumplir que: en X3 se debe cumplir que:
2 12 2 d U 1 (t ) >U 0 (t ) 2 U 1 (t ) U 2 (t )@ ( X 1 ) 2 0,25 2 dt
2 12 2 d U 2 (t ) >U 1 (t ) 2 U 2 (t ) U 3 (t )@ ( X 2 ) 0,25 2 dt 2
2 12 2 d U 3 (t ) >U 2 (t ) 2 U 3 (t ) U 4 (t )@ ( X 3 ) dt 2 0,25 2
Ec.Dif. con Valores de Contorno
0 0
U 4 (t )
en X4 se debe cumplir que: Dr. Ing. A.Mirasso
0
Año 2014
0 13 de 14
Facultad de Ingeniería, Universidad Nacional de Cuyo
Se busca u(x,t) solución de 2 w 2 u ( x, t ) 2 w u ( x, t ) 12 x wx 2 wt 2 u (0, t ) 0
0 en :
^x H R :0 d x d 1`
u(x,t) U1(t)
u ( x,0) sin(S x) wu ( x,0) 3 wt
u (1, t ) 0
U0(t)
0,25
U2(t)
0,5
U3(t) 0,75
U4(t) X=1
Así se puede escribir:
t
2 ª 2 1 0 º U 1 (t ) ½ ª0,25 12 « ° « ° 1 2 1»» ®U 2 (t )¾ « 0 2 « 0,25 «¬ 0 1 2 »¼ °¯U 3 (t ) °¿ «¬ 0
0 0,50 2 0
d 2U 1 (t ) ½ ° ° 2 0 º ° 2dt ° » ° d U 2 (t ) ° 0 »® ¾ dt 2 ° 2 »° 0,75 ¼ d 2U (t ) ° ° 3 ° dt 2 ° ¯ ¿
0½ ° ° ®0¾ °0° ¯ ¿
½ dU 1 ( 0) ° ° U 1 (0) ½ sen(S 0,25)½ ° dt °° °3½° ° ° ° ° ° dU 2 con las condiciones iniciales ®U 2 (0)¾ ®sen(S 0,50) ¾ y ® (0)¾ ®3¾ °U (0) ° °sen(S 0,75) ° ° dt ° °3° ¯ 3 ¿ ¯ ¿ ° dU 3 ¯ ¿ ( 0) ° °¿ °¯ dt Este sistema se puede resolver con métodos analíticos o con métodos numéricos.
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
14 de 14
X
Facultad de Ingeniería, Universidad Nacional de Cuyo
SOLUCIÓN NUMÉRCIA DE ECUACIONES DIFERENCIALES CON VALORES INICIALES
Ecuaciones Diferenciales Ordinarias de Primer Orden Método de Euler Método de Runge Kutta de 2do Orden Sistemas de Ecuaciones Diferenciales Ordinarias de Primer Orden Ecuaciones Diferenciales Ordinarias de Segundo Orden Reducción a Sistemas de primer orden Sistemas de Ecuaciones Diferenciales Ordinarias de Segundo Orden Método de Diferencia Central Reducción a Sistemas de primer orden Ejemplos Resumen
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
1 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
ECUACIONES DIFERENCIALES ORDINARIAS DE PRIMER ORDEN
du (t ) ° f (t , u (t )) ® dt °¯ u (t 0 ) U 0
du (t ) ° A u (t ) ® dt °¯ u (t 0 ) U 0
Se busca obtener una solución aproximada en forma discreta de la siguiente EDO
u(t)
u(t) u0
Forma genércia. Ejemplo La solución discreta es tal que aproxima a la solución exacta del problema
uk # u(tk)
Solución en base a Serie de Taylor en tk se conoce uk y sus derivadas
u1 t0
u (tk 't )
u (tk ) 't
tk
Métodos basados en Derivación Numérica du (t ) 1 (u n1 u n ) O('t ) se consideran aproximaciones dt ' t tn du (t ) dt tn
1 (u n 1 u n 1 ) O('t 2 ) 't
Métodos basados en Integración Numérica
³ du
ub
se considera la integral definida
ua
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
³
t2
t1
du dt
b
f (t , u (t )) dt
u3
u2
' t d 2u 2 dt 2
t
t3
.......... tk
1 (u n 1 u n ) O('t ) 't
du (t ) dt tn 1
ub ua
¦w NP
k
f (tk , uk ) O('t p )
k 1
a
Año 2014
2 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
ECUACIONES DIFERENCIALES ORDINARIAS DE PRIMER ORDEN Se busca obtener una solución aproximada en forma discreta de la siguiente EDO
du (t ) ° f (t , u (t )) ® dt °¯ u (t 0 ) U 0
du (t ) ° A u (t ) ® dt °¯ u (t 0 ) U 0
u(t)
u(t) u0
Solución en base a Serie de Taylor
u1 t0
't 2 d 2u du .......... u (t k 't ) u (t k ) 't dt tk 2 dt 2 t
u2 t2
t1
k
du (t ) dt tk
en tk se conoce uk y
u (t k 't )
f (t k , u (t k ))
f (t k , u k )
fk
d 2 u (t ) dt 2 tk
u3
t
t3
wf (t , u (t )) wf (t , u (t )) du wt wu dt
tk
Solución de Primer Orden
u (t k ) 't f (t k , u k ) O('t 2 )
Solución de Segundo Orden
du 't 2 § wf (t , u (t )) wf (t , u (t )) du · ¸ O('t 3 ) u (t k ) 't ¨ 2 © dt tk wu dt ¹ tk wt
u (t k 't )
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
3 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
EDO DE PRIMER ORDEN-MÉTODO DE EULER Se busca obtener una solución aproximada de la EDO du (t ) ° f (t , u (t )) ® dt °¯ u (t 0 ) U 0 Métodos basados en Derivación Numérica du (t ) 1 (u n1 u n ) O('t ) EULER Adelante considera que dt tn 't
du (t ) dt tn
un1 t n1
un 't f (t n , un ) O('t 2 )
t n 't
tn1 Dr. Ing. A.Mirasso
un+1
f (t n , u n )
u(t) un
du (t ) dt tn 1
1 (u n 1 u n ) O('t ) 't
du (t ) dt tn 1
f (t n 1 , u n 1 )
't fn un
EXPLÍCITO
EULER Atrás considera que
un1
u(t)
t tn
tn+1
un 't f (tn1 , un1 ) O('t 2 )
tn 't
IMPLÍCITO Ec.Dif. con Valores de Contorno
Año 2014
4 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE EULER ATRÁS -EJEMPLO La Ecuación diferencial ordinaria de primer orden
dy (t ) dt
2 t y (t ) 2
u(t)
yex (t ) 1 /(1 t )
u(t)
2
tiene la solución exacta Se busca obtener una solución aproximada de la EDO mediante el Método de EULER ATRÁS Dado (tm, ym),la nueva solución (tm+1, ym+1) se calcula con
y m1
t m1
't fn
un+1
con y (0) 1
y m 't f t m1 , y m1
t m 't
un
un t tn
tn+1
Identificar la función f(t,y). Obtener una solución aproximada en el intervalo >0,1), con el método de Euler Atrás y el paso 't=0.25
…………….se debe resolver una ecuación no lineal!!!!!!!!
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
5 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE EULER ADELANTE -EJEMPLO La Ecuación diferencial ordinaria de primer orden
dy (t ) dt
2 t y (t )
con y (0) 1
yex (t ) 1 /(1 t )
u(t)
tiene la solución
u(t)
2
exacta Se busca obtener una solución aproximada de la EDO mediante el Método de EULER ADELANTE Dado (tm, ym), la nueva solución (tm+1, ym+1) se calcula con
k1 ym1 t m1
't fn
un+1
un
un t
't f t m , ym
tn
tn+1
ym k1
t m 't
Identificar la función f(t,y). Obtener una solución aproximada en el intervalo >0,1@, con el método de Euler adelante y el paso 't=0.25: t y yex(t) Error=abs(yex(t)-yn) k1='t *f(t,y) y(t+'t) 0 1 0,25 0,5 Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
6 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE EULER-EJEMPLO La Ecuación diferencial ordinaria de primer orden
dy (t ) y (t ) dt 2
t 2
con y (0)
4
tiene la solución exacta
Se busca obtener una solución aproximada de la EDO mediante el Método de EULER (ADELANTE) Dado (tm, ym), la nueva solución (tm+1, ym+1) se calcula con
k1 ym1 t m1
yex (t ) 6 e
t / 2
2t
u(t) un+1
't f t m , ym
u(t)
ym k1
un
t m 't
't fn un t
tn+1 tn Identificar la función f(t,y). Obtener una solución aproximada en el intervalo >0,1@, con el método de Euler adelante y el paso 't=0.25: t y yex(t) Error=abs(yex(t)-yn) k1='t *f(t,y) y(t+'t) 0 4 0,25 0,5 0,75 1 Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
7 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
Método de EULER (ADELANTE) en MATLAB
function Euler_00 % Valores Iniciales y0=4; t0=0; Dt=0.01; % incremento de tiempo NDt=1000; % cantidad de Dt a realizar % Dimensionamiento t=zeros(NDt,1); % vector para tiempo y=zeros(NDt,1); % vector para y(t)
dy (t ) y (t ) t dt 2 2 con y (0) 4
Conocido t m , y m
% Inicialización t(1)=t0; y(1)=y0; % EULER for j=1:NDt-1 k1=Dt*(-(1/2)*y(j) + (1/2)*t(j)); y(j+1) = y(j)+k1; t(j+1) = t(j)+Dt; end % Graficación figure(1) plot(t,y(:,1),'b'); end
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
k1 y m 1 t m 1
Año 2014
't f t m , y m y m k1
t m 't
8 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE EULER-EJEMPLO La Ecuación diferencial ordinaria de primer orden
dy (t ) y (t ) dt 2
t) 2
con y (0)
4
tiene la solución exacta
Con 't=0.01 se obtienen las siguientes gráficas
yex (t ) 6 e
t / 2
2t
Solución de EDO -3
9 6
y aprox y exacta
8
x 10
5
7 4
6 3
5 2
4 1
3 0
2
0
Dr. Ing. A.Mirasso
1
2
3
4
5
6
7
8
Ec.Dif. con Valores de Contorno
9
0
1
2
3
4
5
6
7
8
9
10
Año 2014
9 de 20
10
Facultad de Ingeniería, Universidad Nacional de Cuyo
EDO DE PRIMER ORDEN-MÉTODO DE PUNTO MEDIO Se busca obtener una solución aproximada de la EDO du (t ) ° f (t , u (t )) ® dt °¯ u (t 0 ) U 0 Métodos basados en Derivación Numérica Punto Medio
considera que
du (t ) dt tn
1 (u n 1 u n 1 ) O('t 2 ) 2't
du (t ) dt tn
f (t n , u n )
u(t)
u(t)
un+1
un1 tn1
un1 2't f (tn , un ) O('t )
t n 't
3
un-1 un-1
Método Multipaso Orden del Error es 3
Dr. Ing. A.Mirasso
2't fn
un
tn-1
Ec.Dif. con Valores de Contorno
Año 2014
tn
t
tn+1
10 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
EDO DE PRIMER ORDEN-MÉTODOS DE RUNGE KUTTA DE SEGUNDO ORDEN Dado (tm, ym) se busca calcular la nueva solución (tm+1, ym+1)
y m 't ) (t m , y m , 't )
y m1
t m 't
t m1
fm
y(x) G LG
)(tm , ym , 't) a1 f (tm , ym ) a2 f (tG , yG ) tG (tm b1 't)
ym
tm
b1't
Se expande en Taylor
wf wt
tm
(b2't f m )
wf O('t 2 ) wy t
° wf 't )(tm , ym , 't) 't ®(a1 a2 ) f (tm , ym ) a2 (b1't) wt °¯ Dr. Ing. A.Mirasso
m
Ec.Dif. con Valores de Contorno
tG
b2 't
fG b2 't fm
't )
)
A
yG ( ym b2't fm )
(b1't)
Y=Y(x)
ym+1
f (tm , ym )
)
yG
Siendo
f (tG , yG )
L1 E
tm+1
t
't
½ wf 2 ° a2 (b2 't f m ) O('t )¾ wy t tm °¿ m Año 2014
11 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
EDO DE PRIMER ORDEN-MÉTODOS DE RUNGE KUTTA DE SEGUNDO ORDEN Dado (tm, ym) se busca calcular la nueva solución (tm+1, ym+1)
y m1
t m1
y m 't ) (t m , y m , 't )
t m 't
fm
y(x) yG
G LG
)(tm , ym , 't) a1 f (tm , ym ) a2 f (tG , yG ) tG (tm b1 't)
ym
yG ( ym b2't fm )
)
Y=Y(x)
ym+1
Siendo
L1 E
't
)
A
tm
b1't
tm+1
tG b2 't
b2 't
t
't
° ½ wf wf 2 ° 't )(tm , ym , 't) 't ®(a1 a2 ) f (tm , ym ) a2 (b1't) a2 (b2't fm ) O('t )¾ t y w w tm °¯ °¿ tm wf wf a2 (b2 't 2 f m ) O ( 't 3 ) y m1 y m 't (a1 a2 ) f (t m , y m ) a2 (b1't 2 ) wt tm wy t m
't 2 § wf (t , u (t )) wf (t , u (t )) du · du ¸ O('t 3 ) ¨ u (t k 't ) u (t k ) 't wu wt dt tk 2 © dt ¹ tk
Se compara con la solución por Serie de Taylor
Para que las dos Series sean iguales los coeficientes deben ser iguales Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
fG
12 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
EDO DE PRIMER ORDEN-MÉTODOS DE RUNGE KUTTA DE SEGUNDO ORDEN De comparar las series se obtiene que a2 Z z 0 a a 1
a2 b1 1 / 2 a2 b2 1 / 2 1
2
o bien
a1 1 Z
b1
b2
G
t m1
LG
1 2Z
k2 ym1
ym
t m 't
A
tG
't f tG , yG
yG
ym 1 Z k1 Z k2
b1't
t m 't /( 2Z ) 1 ym k1 2Z
tG
b2 't
fG b2 't fm
't )
)
tm
't f xm , y m
t m1
Y=Y(x)
ym+1
y m 't ) (t m , y m , 't )
mediante
k1
)
yG
Dado (tm, ym)se calcula la nueva solución (tm+1, ym+1)
y m1
fm
y(x)
L1 E
tm+1
t
't
t m 't
ª§ 't 't · § ·º ½ y m 't ®1 Z f t m , y m Z f «¨ t m f t m , y m ¸» ¾ O('t 3 ) ¸ , ¨ ym 2Z ¹ © 2Z ¹¼ ¿ ¬© ¯
que es equivalente ha obtener:
y m1 Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
13 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
EDO DE PRIMER ORDEN-MÉTODOS DE RUNGE KUTTA DE SEGUNDO ORDEN Dado (tm, ym), la nueva solución (tm+1, ym+1) se calcula con 't f xm , y m
k1
t m 't /( 2Z )
tG
G
yG
1 yG y m k1 2Z k 2 't f tG , yG t m 't
t m1
tm
b1't
Método de EULER MODIFICADO con Z=1 Dado (tm, ym), la nueva solución (tm+1, ym+1) se calcula con k1 't f t m , y m
tG
y m 1 t m1
Y(t)
y m k1 / 2
't f t G , yG
yG k2
t m 't / 2
L1
t m 't
't
E LC
h)
A
't /2
Ec.Dif. con Valores de Contorno
t
Y=Y(t)
B
tm
Dr. Ing. A.Mirasso
b2 't
tm+1
ll LC
Ym+1 ym+1 Ym
tG
G
YG
ym k 2
't
)
A
ym
)
Año 2014
t tG
't
fG
Y=Y(x ) b2 't fm
LG
ym+1
y m 1 Z k1 Z k 2
y m1
fm
y(x)
L1 E
tm+1
14 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE EULER MODIFICADO Dado (tm, ym), la nueva solución (tm+1, ym+1) se calcula con Z=1 k1 't f t m , y m
tG
y m 1
y m k1 / 2
't f t G , yG
yG k2
t m 't / 2
y(t)
LC
t 2
con y (0)
4
y=y(t)
't ) B
A
t m1 t m 't EJEMPLO La Ecuación diferencial ordinaria de primer orden
dy (t ) y (t ) dt 2
G
yG yexm+1 ym+1 ym
ym k 2
E
L1
tm
tiene la solución exacta
't/2
tG
't
tm+1
ll LC
t
yex (t ) 6 e t / 2 2 t
Se busca obtener una solución aproximada de la EDOmediante el Método de EULER MODIFICADO Obtener una solución aproximada en el intervalo >0,1@, con el método de Euler adelante y el paso 't=0.25: t y k1= tg= yg= k2= Y(t+'t)= 't *f(t,y) t + 't /(2w) y+k1/(2w) 't *f(tG,yG) y+(1-w)k1+w*k2
0 4 0,25 0,5 Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
15 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
Método de EULER MODIFICADO en MATLAB function Euler_00 % Valores Iniciales y0=4; t0=0; Dt=0.01; % incremento de tiempo NDt=1000; % cantidad de Dt a realizar % Dimensionamiento t=zeros(NDt,1); % vector para tiempo y=zeros(NDt,1); % vector para y(t) % Inicialización t(1)=t0; y(1)=y0; w=1
dy (t ) y (t ) t dt 2 2 con y (0) 4
't f t m , y m
k1
% EULER for j=1:NDt-1 k1=Dt*(-(1/2)*y(j) + (1/2)*t(j)); yg = y(j) + k1/(2*w); tg = t(j) + Dt/(2*w); k2=Dt*(-(1/2)*yg + (1/2)*tg); y(j+1) = y(j)+(1-w)*k1+w*k2; t(j+1) = t(j)+Dt; end
tG
y m1 t m1
y m k1 / 2
't f t G , yG
yG k2
t m 't / 2
ym k 2
t m 't
% Graficación figure(1) plot(t,y(:,1),'b'); end Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
16 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE EULER-MODIFICADO La Ecuación diferencial ordinaria de primer orden
dy (t ) y (t ) dt 2
t) 2
con y (0)
4
tiene la solución exacta
Con 't=0.01 se obtienen las siguientes gráficas
yex (t ) 6 e
t / 2
2t
Solución de EDO -5
9 1
y aprox y exacta
8
x 10
0.9 0.8
7
0.7 0.6
6
0.5
5
0.4 0.3
4 0.2 0.1
3
0
2
0
1
Dr. Ing. A.Mirasso
2
3
4
5
6
7
8
9
Ec.Dif. con Valores de Contorno
0
1
2
3
4
5
6
7
8
9
10
10
Año 2014
17 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
EDO DE PRIMER ORDEN-MÉTODOS DE RUNGE KUTTA DE SEGUNDO ORDEN Dado (tm, ym), la nueva solución (tm+1, ym+1) se calcula con 't f xm , y m
k1
t m 't /( 2Z )
tG
k2 y m1
G
yG
LG
ym+1
ym
1 k1 2Z 't f tG , yG
yG
fm
y(x)
y m 1 Z k1 Z k 2
t m 't
t m1
tm
b1't
Método de EULER MEJORADO con Z=1/2 Dado (tm, ym), la nueva solución (tm+1, ym+1) se calcula con k1 't f xm , y m Y(x)
tG
y m 1 t m1
y m k1
't f t G , yG
yG k2
t m 't
y m (k1 k 2 ) / 2
t m 't
L1
Ym
b2 't
tm+1
t
't
ll Lp Y=Y(x)
B A
h) X
Xm Dr. Ing. A.Mirasso
't
L2
E
Ym+h Y’m Ym+1 ym+1
Lp
tG
)
Y=Y(x ) b
)
A
ym
L1 E
Ec.Dif. con Valores de Contorno
h
Año 2014
Xm+1 18 de 20
fG 2 't fm
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE EULER MEJORADO Dado (tm, ym), la nueva solución (tm+1, ym+1) se calcula con Z=1/2 k1 't f xm , y m
tG
y m 1
y m k1
L1
ym+k1
y m (k1 k 2 ) / 2
B A
t m1 t m 't EJEMPLO la EDO
dy (t ) y (t ) dt 2
t 2
tm
con y (0)
4
ll Lp
y=y(t)
ym+1
ym
Lp
L2
E
yexm+1
't f t G , yG
yG k2
t m 't
y(t)
tiene la solución exacta
't
tm+1
't ) t
yex (t ) 6 e t / 2 2 t
Se busca obtener una solución aproximada de la EDOmediante el Método de EULER MEJORADO Obtener una solución aproximada en el intervalo >0,1@, con el método de Euler adelante y el paso 't=0.25: t y k1= xg= yg= k2= Y(t+'t)= 't *f(t,y) t + 't /(2w) y+k1/(2w) 't *f(tG,yG) y+(1-w)k1+w*k2
0 4 0,25 0,5 Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
19 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
EDO DE PRIMER ORDEN-MÉTODO DE LOS TRAPECIOS
du (t ) ° f (t , u (t )) ® dt °¯ u (t 0 ) U 0
Se busca obtener una solución aproximada en forma discreta de la siguiente EDO
u(t)
u(t) un-1
Métodos basados en Integración Numérica se considera la integral definida
³ du
un 1
un
un1
³
tn 1
f (t , u (t )) dt
un tn-1
tn
un+1
t
tn+1
't un f (tn , un ) f (tn , un1 ) O('t 3 ) 2 tn
Método Implícito de un paso Solución de ecuación no lineal Corrector iterativo de una Predicción con algún método explícito
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
20 de 20
Facultad de Ingeniería, Universidad Nacional de Cuyo
SOLUCIÓN NUMÉRCIA DE ECUACIONES DIFERENCIALES CON VALORES INICIALES
Ecuaciones Diferenciales Ordinarias de Primer Orden Método de Euler Método de Runge Kutta de 2do Orden Sistemas de Ecuaciones Diferenciales Ordinarias de Primer Orden Ecuaciones Diferenciales Ordinarias de Segundo Orden Reducción a Sistemas de primer orden Sistemas de Ecuaciones Diferenciales Ordinarias de Segundo Orden Método de Diferencia Central Reducción a Sistemas de primer orden Ejemplos Resumen
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
1 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
SISTEMAS DE EDO -CIRCUITOS ELÉCTRICOS Circuito Butterworth para filtros pasa baja (Eronini Umez Eronini, Dinàmica de Sistemas y Control, Thomson Learning, 2001). Se trata de un circuito de tipo RLC; es decir, con resistencias R, VB VC VD inductancias L, y capacitares C. Las ecuaciones que describen el comportamiento son las siguientes: RS L2 C2 di1 ½ L1 ° dt ° VD(t) i1(t) RL Ve(t) i2(t) 0 0 º i1 ½ Ve / L1 ½ ° dV ° ª RS / L1 1 / L1 C 1 » °V ° ° 0 ° ° C ° « 1/ C 0 0 1/ C2 ° dt ° « ° ° C° ° 1 » ® di ¾ ¾ ® ¾® « » 0 1 / 0 1 / L L 2 2 2 ° ° ° i2 ° ° 0 ° » ° dt ° «¬ 0 0 1 / C 2 1 / RL C 2 ¼ °¯VD °¿ °¯ 0 °¿ ° dVD ° °¯ dt °¿ Circuitos RL y RC (Zill,G. Cullen, M. Ecuaciones Diferenciales con Problemas de Valores de Frontera, Thomson Learning, 2002). Se trata de circuitos de tipo RL; es decir, con resistencias R, inductancias L, y RC con resistencias y capacitores C. Cuando se consideran como incógnitas las corrientes, se tiene que las ecuaciones resultantes son las siguientes: i2(t)
L1
di2 (t ) ( R1 R2 ) i2 (t ) R1 i3 (t ) Ve(t ) dt di (t ) L2 3 R1 i2 (t ) R1 i3 (t ) Ve(t ) dt i1 (t ) i2 (t ) i3 (t )
R1
Ec.Dif. con Valores de Contorno
i3(t) L1 R2
Ve(t)
di (t ) L 1 R i2 (t ) Ve(t ) dt di (t ) RC 2 i2 (t ) i1 (t ) 0 dt i1 (t ) i2 (t ) i3 (t ) Dr. Ing. A.Mirasso
i1(t)
i1(t) L R
L2
i3(t) i2(t) C
Ve(t)
Año 2014
2 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
SISTEMAS DE EDO-MODELOS DE DISPERSIÓN Modelo de Dispersión de antihistaminas. (Borelli, R.; Colleman, C. Ecuaciones Diferenciales, Oxford University Press, 2002). Durante los resfríos es típico ingeniar sustancias para aliviar los malestares. Dichas sustancias, como la antihistaminas, se presentan en el organismo como sustancias a eliminar. En principio las antihistaminas están en el aparato digestivo, y desde allí pasan a la sangre que las circula por todo el cuerpo y son eliminadas de la sangre en los riñones. Se considera con x(t) es la cantidad de antihistamina en el aparato digestivo en el instante de tiempo t; y con y(t) se considera la cantidad de antihistamina en la sangre en el mismo instante de tiempo t. Así un modelo de intercambio de antihistaminas en el organismo está dado por las siguientes dos consideraciones: x La tasa de cantidad de antihistaminas en el aparato digestivo es proporcional a la cantidad de antihistaminas en el aparato digestivo. Con una cierta cantidad inicial de antihistaminas
dx (t ) k1 x (t ) dt x ( 0) x 0
k1 x(t) x(t)
k2 y(t) y(t)
Al ser una tasa negativa, deja en manifiesto que el proceso es tal que la cantidad decrece. x La tasa de cantidad de antihistaminas en sangre es proporcional a la cantidad de antihistaminas en sangre y a la cantidad de antihistaminas que ingresan desde el aparato digestivo. Con una cierta cantidad inicial de antihistaminas dy (t ) k 2 y (t ) k1 x (t ) dt y ( 0) y 0 Otra forma de expresar las ecuaciones diferenciales es en la siguiente forma matricial. dx(t ) ½ ° dt ° ® dy (t ) ¾ ° ° ¯ dt ¿
Dr. Ing. A.Mirasso
ª k1 «k ¬ 1
0 º x(t ) ½ ¾ ® k 2 »¼ ¯ y (t )¿
Ec.Dif. con Valores de Contorno
con
x ( 0) ½ ¾ ® ¯ y ( 0) ¿
Año 2014
x0 ½ ® ¾ ¯ y0 ¿
3 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo Eliminación de contaminante en un líquido Considérese un proceso por el cual se busca eliminar un contaminante de un liquido que circula entre dos depósitos depuradores. Se trata de dos bloques o sistemas depuradores en los cuales la cantidad de contamínate en el primer bloque es x(t), mientras que en el segundo es y(t). En la siguiente figura se muestra el diagrama de bloques del sistema en análisis. En el primer bloque solo salen k1 x(t) cantidad de contaminante, que ingresan en el k1 x(t) k2 y(t) segundo bloque. Pero también existe una retroalimentación no deseada desde el segundo x(t) y(t) bloque que hace ingresar al primer bloque k3 y(t) contaminante. Desde el segundo bloque k3 y(t) salen k2 y(t) cantidad de contaminantes residuales. Las cantidades que salen de un bloque se consideran negativas, mientras que las que ingresan, positivas. De esa manera es posible hacer el balance de lo que entra y sale en cada bloque e igualarlo a la tasa de cantidad de sustancia contaminante en el tiempo de dicho bloque.
dx(t ) k 3 y (t ) k1 x(t ) dt x ( 0) x 0 dx(t ) ½ ° dt ° ® dy (t ) ¾ ° ° ¯ dt ¿
dy (t ) k 3 y (t ) k1 x(t ) k 2 y (t ) dt y (0) y 0
Otra forma de expresar las ecuaciones diferenciales es en la siguiente forma matricial.
ª k1 «k ¬ 1
º x(t ) ½ ® ¾ k 2 k 3 »¼ ¯ y (t )¿ k3
con
x(0) ½ ¾ ® ¯ y (0)¿
x0 ½ ® ¾ ¯ y0 ¿
En este caso se debe resolver las dos incógnitas en forma simultánea con las dos ecuaciones.
Otros problemas descriptos por SISTEMAS de EDO son: Flujos entre depósitos; Flujos en procesos químicos; Evolución de poblaciones y/o especies Evolución de “Romeo y Julieta”
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
4 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
SISTEMAS DE ECUACIONES DIFERENCIALES ORDINARIAS DE PRIMER ORDEN Se busca obtener una solución aproximada en forma discreta del siguiente SISTEMA de EDO
dy1 dt dy 2 dt
10 y1 (t ) 4 y 2 (t ) 4 y1 (t ) 0 y 2 (t )
, con:
y1 (0)
5
y 2 ( 0)
3 , que tiene solución exacta
Aplicar el Método de Euler y obtener la siguiente aproximación t 0,01 0,02 0,03 0 y1 5 y2 3 k1 _y1 0,01*(-10*5+4*3) _y2 0,01*( -4*5+0*3) y1_(n+1) 5+(-0,38) y2_(n+1) 3+(-0,2) tn+1 0+0,1
y (t )
1 1½ 2*t 14 1 ½ 8*t ® ¾e ® ¾e 3 ¯1 / 2¿ 3 ¯2¿
u(t) 0,04
0,05
un+1
't fn
u(t) un
un t tn+1
tn
't f t m , y m
Método de Euler: Dado (tm,ym)
k1 y m 1 t m 1 Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
y m k1
t m 't
5 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
SISTEMAS DE ECUACIONES DIFERENCIALES ORDINARIAS DE PRIMER ORDEN Se busca obtener una solución aproximada en forma discreta del siguiente SISTEMA de EDO
dy1 dt dy 2 dt
10 y1 (t ) 4 y 2 (t ) 4 y1 (t ) 0 y 2 (t )
, con:
y1 (0)
5
y 2 ( 0)
3,
Aplicar el Método de Euler Modificado (w=1) y obtener la siguiente aproximación t 0 0,01 0,02 0,03 0,04 y1 5 4,635 4,2977 y2 3 2,8076 2,6292 k1 _y1 0,01*(-10*5+4*3)=-0,38 _y2 0,01*( -4*5+0*3)=-0,2 tg yg1 yg2 k2 _y1 _y2
't f t m , y m
k1
tG
y m 1
0+0,01/2=0,005 5 -0,38/2= 4,81 3 -0,2/2= 2,9 0,01*(-10*4,81+4*2,9)= -0,365 0,01*( -4*4,81+0*2,9)= -0,1924
t m 1
ym
1 k1 2Z 't f tG , yG
yG k2
t m 't /( 2Z )
y m 1 Z k1 Z k 2
t m 't
tg 0+0,01 y1_(n+1) 5 -0,365= 4,635 y2_(n+1) 3 -0,1924= 2,8076 Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
6 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo % Datos y10=5; % Valores Iniciales y20=3; t0=0; Dt=0.01; % incremento de tiempo NDt=200; % cantidad de Dt a realizar % Dimensionamiento t=zeros(1,NDt); % vector fila para el tiempo y=zeros(2,NDt); % Matriz para el vector solución y(t) yg=zeros(2,1); % Vector columna auxiliar ya=zeros(2,1); k1=zeros(2,1); k2=zeros(2,1); % Inicialización del primer estado solución t(1)=t0; y(1,1)=y10; y(2,1)=y20; % Euler % Euler Modificado for j=1:NDt-1 for j=1:NDt-1 ya=y(:,j); ya=y(:,j); ta=t(j); ta=t(j); k1(1,1)= Dt* (-10*ya(1) + 4*ya(2)); k1(2,1)= Dt* (-4*ya(1)+0*ya(2)); y(:,j+1) t(j+1) end
= ya + k1; = ta + Dt;
figure(1) plot(t,y(1,:),'b');grid on end
k1(1,1)= Dt* (-10*ya(1) + 4*ya(2)); k1(2,1)= Dt* (-4*ya(1)+0*ya(2)); yg tg
k2(1,1)= Dt* (-10*yg(1) + 4*yg(2)); k2(2,1)= Dt* (-4*yg(1)+0*yg(2)); y(:,j+1) t(j+1) end
for i=1:NDt yex(i)=(1/3)*exp(-2*t(i))+(14/3)*exp(-8*t(i)); end normer=norm(er,inf) figure(1) plot(t,y(1,:),'b', t,yex(1,:),'r');grid on end
Dr. Ing. A.Mirasso
= ya + k1/(2); = ta + Dt/(2);
Ec.Dif. con Valores de Contorno
= ya + k2; = ta + Dt;
er(i)=abs(yex(i)-y(1,i));
Año 2014
7 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo % Datos y10=5; % Valores Iniciales y20=3; t0=0; Dt=0.01; % incremento de tiempo NDt=200; % cantidad de Dt a realizar % Dimensionamiento t=zeros(1,NDt); % vector fila para el tiempo y=zeros(2,NDt); % Matriz para el vector solución y(t) yg=zeros(2,1); % Vector columna auxiliar ya=zeros(2,1); k1=zeros(2,1); k2=zeros(2,1); % Inicialización del primer estado solución t(1)=t0; y(1,1)=y10; y(2,1)=y20; % Euler % Euler Modificado for j=1:NDt-1 for j=1:NDt-1 ya=y(:,j); ya=y(:,j); ta=t(j); ta=t(j); k1=Dt*f_sist_1(ya,ta) k1=Dt*f_sist_1(ya,ta) k1(1,1)= Dt* (-10*ya(1) + 4*ya(2)); k1(1,1)= Dt* (-10*ya(1) + 4*ya(2)); k1(2,1)= Dt* (-4*ya(1)+0*ya(2)); k1(2,1)= Dt* (-4*ya(1)+0*ya(2)); y(:,j+1) = ya + k1; yg = ya + k1/(2); t(j+1) = ta + Dt; tg = ta + Dt/(2); end k2=Dt*f_sist_1(yg,tg); k2(1,1)= Dt* (-10*yg(1) + 4*yg(2)); figure(1) k2(2,1)= Dt* (-4*yg(1)+0*yg(2)); plot(t,y(1,:),'b');grid on y(:,j+1) = ya + k2; t(j+1) = ta + Dt; end end end function [fy]=f_sist_1(z,x) fy(1,1)=(-10*z(1) + 4*z(2)); fy(2,1)=(-4*z(1)+0*z(2)); end
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
% Euler y RK for j=1:NDt-1 ya=y(:,j); ta=t(j); k1=Dt*f_sist_1(ya,ta) if (w==0) k2=zeros(2,1) else yg = ya + k1/(2*w); tg = ta + Dt/(2*w); k2=Dt*f_sist_1(yg,tg); end y(:,j+1) = ya + (1-w)*k1+w* k2; t(j+1) = ta + Dt; end
8 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
Se busca obtener una solución aproximada en forma discreta del siguiente SISTEMA de EDO
10 y1 (t ) 4 y 2 (t )
dy1 dt dy 2 dt
4 y1 (t ) 0 y 2 (t )
y1(t) con Dt=0.02
, con:
y1 (0)
5
y 2 ( 0)
3 , que tiene solución exacta
EULER MODIFICADO
y (t )
1 1½ 2*t 14 1 ½ 8*t ® ¾e ® ¾e 3 ¯2¿ 3 ¯1 / 2¿
EULER
Solución de EDO 5
Solución de EDO 5
y aprox y exacta
4.5
y aprox y exacta
4.5
4 4
3.5 3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
0
2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Soluciones Aproximadas -3
9
x 10
0.16
8
0.14
7
0.12
6
0.1
5
0.08 4
0.06 3
0.04
2
0.02
1 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Funciones Error Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
9 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
ECUACIONES DIFERENCIALES ORDINARIAS DE SEGUNDO ORDEN REDUCCIÓN A SISTEMAS DE PRIMER ORDEN Considérese la EDO de segundo orden del péndulo simple, d 2T (t ) g T (t ) 0 2 L dt
T (t0 ) T 0
dT (t ) dt
E 0 en t t0
donde T(t) es la posición angular medida respecto de la vertical; T0 la posición inicial y E0 la velocidad inicial. Si se plantea un cambio de variables tal que y1 (t ) T (t ) y1 (t ) y2 (t )
y2 (t )
dT (t ) dt
entonces y (t ) 2
Se tiene y 1 (t ) y 2 (t )
y 2 (t )
( g / L) y1 (t )
con lo que
g ( ) T (t ) 0 y t con lo que la EDO es 2 L
y1 (t ) 0 y1 (t ) 1 y2 (t ) y 2 (t ) ( g / L) y1 (t ) 0 y2 (t )
T (t )½ y1 (t ) ½ ® ¾ ® ¾ y t ( ) t ( ) T ¿ ¯ 2 ¿ ¯ & & y (t ) ½ 0 y1 (t ) 1 y 2 (t ) ½ T (0)½ y1 (0) ½ T 0 ½ & f (t , y ) ® 1 ¾ ® ¾ con y (0) ® ¾ ® ¾ ® ¾ ¯ y 2 (t )¿ ¯ ( g / L) y1 (t ) 0 y 2 (t )¿ ¯T (0)¿ ¯ y 2 (0)¿ ¯E 0 ¿
& y (t )
y Dr. Ing. A.Mirasso
o bien
d 2T (t ) dt 2
Ec.Dif. con Valores de Contorno
Año 2014
10 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
En un sistema masa resorte los desplazamientos son la solución de la siguiente EDO de segundo orden: Buscar u(t), para t H >0;f , tal que u (0)½ 0½ d 2 u (t ) du (t ) m c k u t p sen t con : ( ) ( ) ¾ ® ¾ ® dt dt 2 ¯u (0)¿ ¯0¿ Si se plantea un cambio de variables tal que y1 (t ) u (t ) dy1 (t ) du (t ) lo que implica que y 2 (t ) y 2 (t ) dt dt
y la EDO se puede escribir m
dy 2 (t ) c y 2 (t ) k y1 (t ) dt
p sen(: t )
con
El problema original se puede reemplazar por Buscar y1 (t ) , y2 (t ) , para t H >0;f , tal que dy1 (t ) dt dy 2 (t ) dt
y 2 (t )
con y1 (0)
( p / m) sen(: t ) (k / m) u (t )
y 2 (0)
0 0
d u (t )½ ¾ ® dt ¯v(t ) ¿
(: / Z ) sen(Z t ) sen(: t )½ 2 ® ¾ con Z ¯ : ( cos(Z t ) cos(: t )) ¿
La solución exacta de este problema es u (t )½ ¾ ® ¯v(t ) ¿
Dr. Ing. A.Mirasso
p 1 k 1 (: / Z ) 2
Ec.Dif. con Valores de Contorno
k/m
Año 2014
y 1 ( 0) ½ ® ¾ ¯ y 2 (t ) ¿
0½ ® ¾ ¯0¿
v(t ) ½ ¾ con ® ¯( p / m) sen(: t ) (k / m) u (t )¿
u (0)½ ¾ ® ¯v(0) ¿
0½ ® ¾ ¯0¿
t H >0;f
11 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
En un sistema la incógnita primaria es la solución de la siguiente EDO de tercer orden: Buscar u(t), para t H >0;f , tal que du (t ) d 2 u (t ) d 3u (t ) a a0 u (t ) b(t ) a a3 1 2 dt dt 2 dt 3 Si se plantea un cambio de variables tal que y1 (t ) u (t )
u (0)½ ° ° con ®u (0)¾ °u(0)° ¿ ¯
A½ ° ° ®B¾ °C ° ¯ ¿
dy1 (t ) dt dy 2 (t ) dt
y 2 (t ) du (t ) y 2 (t ) lo que implica que dt y3 (t ) d 2 u (t ) y3 (t ) dt 2 dy (t ) a3 3 a2 y3 (t ) a1 y 2 (t ) a0 y1 (t ) b(t ) y la EDO se puede escribir dt
con u (0)
A; u (0)
B; u(0)
C
y2 (t ) , y3 (t ) para t H >0;f , tal que
El problema original se puede reemplazar por:
Buscar y 1 (t ) , dy1 (t ) dt dy 2 (t ) dt dy3 (t ) dt Dr. Ing. A.Mirasso
y 2 (t ) y3 (t ) (a 2 / a3 ) y3 (t ) (a1 / a3 ) y 2 (t ) (a0 / a3 ) y1 (t ) b(t ) / a3 Ec.Dif. con Valores de Contorno
Año 2014
y1 (0) ½ ° ° con ® y 2 (0) ¾ ° y ( 0) ° ¯ 3 ¿
A½ ° ° ®B¾ °C ° ¯ ¿ 12 de 13
Buscar y 1 (t ) , y2 (t ) , y3 (t ) para t H >0;f , tal que dy1 (t ) ½ ° dt ° 0 1 0 º y1 (t ) ½ 0 °° dy (t ) °° ª « » ° y (t )° ° 0 2 0 0 1 ¾ « ® »® 2 ¾ ® dt ° « (a / a ) (a / a ) (a / a )» ° y (t ) ° °b(t ) / a ° dy t ( ) 0 3 1 3 2 3 ¼¯ 3 3 ¿ ¯ ° 3 ° ¬ °¯ dt °¿
Facultad de Ingeniería, Universidad Nacional de Cuyo
½ ° ¾ ° ¿
y1 (0) ½ ° ° con ® y 2 (0)¾ ° y ( 0) ° ¯ 3 ¿
A½ ° ° ®B¾ °C ° ¯ ¿
Ejemplo: Pag. 149 de Zill,G. Cullen, M. Ecuaciones Diferenciales con Problemas de Valores de Frontera, Thomson Learning, 2002. Buscar u(t) solución de
d 3u (t ) d 2 u (t ) du(t ) 6 11 6 u (t ) 3 t dt dt 3 dt 2
con u (0)
A ; u (0)
B u(0) C
c1 e t c2 e 2t c3 e 3t (11 / 2) (1 / 2) t
La solución exacta es:
uex (t )
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
13 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
SOLUCIÓN NUMÉRCIA DE ECUACIONES DIFERENCIALES CON VALORES INICIALES
Ecuaciones Diferenciales Ordinarias de Primer Orden Método de Euler Método de Runge Kutta de 2do Orden Sistemas de Ecuaciones Diferenciales Ordinarias de Primer Orden Ecuaciones Diferenciales Ordinarias de Segundo Orden Reducción a Sistemas de primer orden Sistemas de Ecuaciones Diferenciales Ordinarias de Segundo Orden Método de Diferencia Central Reducción a Sistemas de primer orden Ejemplos Resumen
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
1 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
ECUACIONES DIFERENCIALES ORDINARIAS DE SEGUNDO ORDEN REDUCCIÓN A SISTEMAS DE PRIMER ORDEN En un sistema la incógnita primaria es la solución de la siguiente EDO de segundo orden: Buscar u(t), para t H >0;f , tal que
d 2 u (t ) du(t ) 12 6 48 u (t ) 36 sen(5 t ) 2 dt dt u (t )½ ® ¾ v t ( ) ¯ ¿
(5 / 2) sen(2 t ) sen(5 t )½ 36 1 ® ¾ 48 1 (5 / 2) 2 ¯ 5 ( cos(2 t ) cos(5 t )) ¿
La solución exacta de este problema es
t H >0;f
u (0)½ 0½ con ® ¾ ® ¾ ( 0 ) u ¯ ¿ ¯0¿
Reducir a un sistema de EDO de primer orden de la forma:
& d y (t ) dt
® ¯
expresar
& y (t )
½ ¾ ¿
& y ( 0)
& y0
& & f ( y (t ), t )
® ¯
½ ¾ ¿
& con y (0)
& & f ( y (t ), t )
& y0
® ¯
½ ¾ ¿
y resolver con un método de Runge Kutta de segundo orden. Comparar con la solución exacta. Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
2 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
REDUCCIÓN A SISTEMAS DE PRIMER ORDEN En un sistema no lineal la incógnita primaria es la solución de la siguiente EDO de segundo orden: Buscar u(t), para t H >0;f , tal que
d 2 u (t ) du(t ) c k (1 a u (t )) u (t ) m 2 dt dt
u (0)½ 0½ con ® ¾ ® ¾ ¯u (0)¿ ¯0¿
p sen(: t )
Reducir a un sistema de EDO de primer orden de la forma:
& dy (t ) dt
® ¯
expresar
& y (t )
½ ¾ ¿
& y ( 0)
& y0
& & f ( y (t ), t )
® ¯
½ ¾ ¿
& con y (0)
& & f ( y (t ), t )
& y0
® ¯
½ ¾ ¿
y resolver con un método de Runge Kutta de segundo orden. Comparar con la solución exacta.
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
3 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
REDUCCIÓN A SISTEMAS DE PRIMER ORDEN En un sistema la incógnita primaria es la solución de la siguiente EDO de tercer orden: Buscar u(t), para t H >0;f , tal que d 3u (t ) d 2 u (t ) du(t ) con u (0) A ; u (0) B u(0) C 6 11 6 u (t ) 3 t dt dt 2 dt 3 La solución exacta de este problema es
uex (t )
c1 e t c2 e 2t c3 e 3t (11 / 2) (1 / 2) t
Ejemplo: Pag. 149 de Zill,G. Cullen, M. Ecuaciones Diferenciales con Problemas de Valores de Frontera, Thomson Learning, 2002. Buscar u(t) solución de
Reducir a un sistema de EDO de primer orden de la forma:
& dy (t ) dt
° ® ° ¯
expresar
& y (t )
½ ° ¾ ° ¿
& y ( 0)
& y0
& & f ( y (t ), t )
° ® ° ¯
½ ° ¾ ° ¿
& con y (0)
& & f ( y (t ), t )
° ® ° ¯
& y0
½ ° ¾ ° y ¿
resolver con un método de Runge Kutta de segundo orden, eligiendo Dt y el coeficiente w. Comparar con la solución exacta.
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
4 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
REDUCCIÓN A SISTEMAS DE PRIMER ORDEN En un sistema las incógnitas primarias son la solución del siguiente sistema de EDO de segundo orden: Buscar T(t), T(t) para t H >0;f , tal que
d 2T1(t) 16 4(T 2 T1 ) 5t dt 2 T1 ½ ® ¾ d 2T 2(t) con ¯T 2 ¿ 2 t0 32 9T1 3t 2 dt
dT 1 ½ S / 8 ½ ° dt ° ¾ y ® dT ¾ ® S / 4 ¿ ° 2° ¯ ¯ dt ¿ t 0 Reducir a un sistema de EDO de primer orden de la forma: & d y (t ) & & & & f ( y (t ), t ) con y (0) y 0 dt Expresar
& y (t )
& y ( 0)
& y0
0,1 ½ ¾ ® ¯ 0,2¿
& & f ( y (t ), t )
y resolver con un método de Runge Kutta de segundo orden. Comparar con la solución exacta.
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
5 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
REDUCCIÓN A SISTEMAS DE PRIMER ORDEN En un sistema las incógnitas primarias son la solución del siguiente sistema de EDO de segundo orden: Buscar x1(t), x2(t) x3(t) para t H >0;f , tal que
M x (t ) C x (t ) K x (t )
R (t ),
ª1 0 0º x1 (t ) ½ ª 4 0 0º x1 (t ) ½ ª0 4 1º x1 (t ) ½ «0 1 0» ° x (t ) ° «0 1 0» ° x (t ) ° «4 2 0» ° x (t ) ° »® 2 ¾ « »® 2 ¾ « »® 2 ¾ « ° ° ° ° «¬0 0 2»¼ ¯ x3 (t ) ¿ «¬0 0 3»¼ ¯ x 3 (t ) ¿ «¬1 0 1»¼ °¯ x 3 (t ) °¿
x 1 ( 0) ½ ° ° con ® x 2 (0) ¾ ° x ( 0) ° ¯ 3 ¿
1 ½ ° ° ®2 ¾ °1 ° ¯ ¿
x1 (0) ½ ° ° ® x 2 (0) ¾ ° x (0) ° ¯ 3 ¿
1 ½ ° ° ®4 ¾ °0 ° ¯ ¿
5e t 8e 2t cos t ½ ° ° 2t t ® 8e 4 e ¾ ° cos t 3sent e t ° ¯ ¿
Reducir a un sistema de EDO de primer orden de la forma:
& d y (t ) dt
& & f ( y (t ), t )
& con y (0)
& y0
Expresar
& y (t )
& y ( 0)
& y0
& & f ( y (t ), t )
y resolver con un método de Runge Kutta de segundo orden. Comparar con la solución exacta. Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
6 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
SISTEMA DE EDO DE PRIMER ORDEN -MATLAB % Datos Dim= 2 % Dimensión del Sistema de EDO w= 1 % 0 es Euler; 1 es E-Modificado; ½ es E-Mejorado y0=[5; 3]; % Valores Iniciales deben ser Dim t0=0; Dt=0.01; % incremento de tiempo NDt=200; % cantidad de Dt a realizar % Dimensionamiento t=zeros(1,NDt); % vector fila para el tiempo y=zeros(Dim,NDt); % Matriz para el vector solución y(t) yg=zeros(2,1); % Vector columna auxiliar ya=zeros(Dim,1); k1=zeros(Dim,1); k2=zeros(2,1); % Inicialización del primer estado solución t(1)=t0; y(:,1)=y0; % Euler y RK for j=1:NDt-1 ya=y(:,j); ta=t(j); k1=Dt*f_sist_1(ya,ta) if (w==0) k2=zeros(Dim,1) else yg = ya + k1/(2*w); tg = ta + Dt/(2*w); k2=Dt*f_sist_1(yg,tg); end y(:,j+1) = ya + (1-w)*k1+w* k2; t(j+1) = ta + Dt; end end function [fy]=f_sist_1(z,x) % deben ser Dim componentes fy(1,1)=(-10*z(1) + 4*z(2)); fy(2,1)=(-4*z(1)+0*z(2)); end
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
7 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
ECUACIONES DIFERENCIALES ORDINARIAS DE SEGUNDO ORDEN MÉTODO DE DIFERENCIA CENTRAL En un sistema la incógnita primaria es la solución de la siguiente EDO de segundo orden: Buscar u(t), para t H >0;f , tal que
d 2u (t ) du(t ) k u (t ) m c 2 dt dt
u (0)½ con ® ¾ u ( 0 ) ¯ ¿
p g (t )
En base a las reglas de derivadas centrales d 2 u (t ) 1 (u n 1 2 u n u n 1 ) O('t 2 ) 2 2 't dt t
u(t)
1 (u n 1 u n 1 ) O('t 2 ) 2 2 't
un
Se puede escribir la EDO en tn en la forma:
un-1
u n1
bg (t n ) DE u n H E u n1
Con 't · § ¨m c¸ ; 2 ¹ © 1
GE
DE
G E 2 m 't 2 k ;
u(t)
un+1
n
du (t ) dt t n
D ½ ® ¾ ¯E ¿
HE
tn-1
't
· § 't GE ¨ c m ¸; bg (t ) ¹ © 2
tn
't
t tn+1
' t 2 G E p g (t )
Es necesario conocer dos estados solución (un-1; un) para calcular el nuevo estado un+1 en tn+1. Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
8 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
EJEMPLO Buscar con el método de diferencia central la incógnita u(t), para t H >0;f , tal que
d 2u (t ) 65 400 u (t ) 80 sen(5 t ) 2 dt
80 1 ^ (5 / 2.48) sen(2.48 t ) sen(5 t )` 400 1 (5 / 2.48) 2
La solución exacta de este problema es u (t )
u (0)½ con ® ¾ ( 0 ) u ¯ ¿
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
0½ ® ¾ ¯0¿
9 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
MÉTODO DE DIFERENCIA CENTRAL
(t ) C u (t ) K u (t ) R (t ), Mu
Para un sistema EDO que en forma generalizada se puede escribir
con valores iniciales conocidos u (0), u (0), y con u (0) que se puede obtener de satisfacer la ecuación diferencial vectorial en el instante inicial t=0. (t ) y la u (t ) con derivadas numéricas centrales y se reemplazan en el sistema original. Se aproximan la u M
1 u (t 't ) 2u (t ) u (t 't ) C 1 u (t 't ) u (t 't ) Ku (t ) 2 't 2 't
Es posible aproximar u (t 't ) , para cualquier t t t0 en la forma:
u (t ' t )
con
b g (t ) D E ( ' t ) u ( t ) H E ( ' t ) u ( t ' t )
't · § ¨M C¸ ; 2 ¹ © 1
GE
R (t )
DE
G E 2 M 't 2 K ;
HE
§ 't · G E ¨ C M ¸; b g (t ) © 2 ¹
' t 2 G E R (t )
Se debe destacar que las matrices GE y HE, como así también la matriz inversa que las define, se calcula sólo una vez al inicio del método. Para el valor inicial de t (t=t0), primero se halla una aproximación de u (t0 't ) mediante Serie de Taylor: u (t 0 't )
Dr. Ing. A.Mirasso
1 2 § · (t 0 ) ¸ . ¨ u (t 0 ) 't u (t 0 ) 't u 2 © ¹
Ec.Dif. con Valores de Contorno
Año 2014
10 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo
EJEMPLO Buscar con el método de diferencia central las incógnitas x1(t), x2(t) x3(t) para t H >0;f , tal que
M x (t ) C x (t ) K x (t )
R (t ),
ª1 0 0º x1 (t ) ½ ª 4 0 0º x1 (t ) ½ ª0 4 1º x1 (t ) ½ «0 1 0» ° x (t ) ° «0 1 0» ° x (t ) ° «4 2 0» ° x (t ) ° »® 2 ¾ « »® 2 ¾ « »® 2 ¾ « «¬0 0 2»¼ °¯ x3 (t ) °¿ «¬0 0 3»¼ °¯ x 3 (t ) °¿ «¬1 0 1»¼ °¯ x 3 (t ) °¿
x 1 ( 0) ½ ° ° con ® x 2 (0) ¾ ° x ( 0) ° ¯ 3 ¿
Dr. Ing. A.Mirasso
1 ½ ° ° ®2 ¾ °1 ° ¯ ¿
x1 (0) ½ ° ° ® x 2 (0) ¾ ° x (0) ° ¯ 3 ¿
Ec.Dif. con Valores de Contorno
1 ½ ° ° ®4 ¾ °0 ° ¯ ¿
Año 2014
5e t 8e 2t cos t ½ ° ° 2t t ® 8e 4 e ¾ t ° cos t 3sent e ° ¯ ¿
11 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo function Dif_cen clc,clear % Datos Dt=0.003; % incremento de tiempo NDt=2500; % cantidad de Dt a realizar dim=3; % cantidad de funciones incógnitas % Dimensionamiento %t=zeros(1,NDt); %y=zeros(dim,NDt); %yan=zeros(dim,1); %yac=zeros(dim,1); %ynu=zeros(dim,1); M= [1 0 0 C= [4 0 0 K= [0 4 1
0 -1 0 0 -1 0 4 2 0
% % % % %
vector Matriz Vector Vector Vector
fila para guardar el tiempo para guardar el vector solución y(t) de solución anterior de solución actual de solución nueva
0 0 2]; 0 0 3]; 1 0 1];
% Valores Iniciales o actuales tac=0; yac= [1; 2; 1] vac= [1; 4; 0] fua(1,1)=5*exp(tac)+8*exp(2*tac)+cos(tac); fua(2,1)=-8*exp(2*tac)+4*exp(tac); fua(3,1)=-cos(tac)-3*sin(tac)+exp(tac); % Inicialización G=inv(M+(Dt/2)*C) D=G*(2*M-Dt^2*K) H=G*((Dt/2)*C-M) yan=yac-Dt*vac+((Dt^2)/2)*inv(M)*(fua-K*yac-C*vac); t(1)=tac; y(:,1)=yac;
Dr. Ing. A.Mirasso
% Almacenamiento para luego graficar
Ec.Dif. con Valores de Contorno
Año 2014
12 de 13
Facultad de Ingeniería, Universidad Nacional de Cuyo % Diferencia Central for j=2:NDt bac=fun_ind(tac,G,Dt); % actualización de término independiente ynu=bac + D*yac + H*yan; % Calculo con Dif Central tnu=tac+Dt; t(j)=tnu; % Almacenamiento para luego graficar y(:,j)=ynu; yan=yac; yac=ynu; tac=tnu; end end
% actualización de estado anterior % actualización de estado actual
function [fy]=fun_ind(x,G,Dt) r(1,1)= 5*exp(x) +8*exp(2*x) +cos(x); r(2,1)=-8*exp(2*x)+4*exp(x); r(3,1)=-cos(x) -3*sin(x) +exp(x); fy=Dt^2*G*r; end
Dr. Ing. A.Mirasso
Ec.Dif. con Valores de Contorno
Año 2014
13 de 13