Cnyc_2016_material_para_clases_teoria.pdf

  • Uploaded by: Augusto
  • 0
  • 0
  • November 2019
  • PDF

This document was uploaded by user and they confirmed that they have the permission to share it. If you are author or own the copyright of this book, please report to us by using this DMCA report form. Report DMCA


Overview

Download & View Cnyc_2016_material_para_clases_teoria.pdf as PDF for free.

More details

  • Words: 27,737
  • Pages: 164
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 a˜x b˜xc 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

More Documents from "Augusto"