Post on 02-Feb-2016
transcript
/48
Simulacion de sistemas dinamicos
Principios básicos de la integración numérica
1
/48
Contenido
El dominio de estabilidad numérica
Cálculo del dominio de estabilidad numérica
Iteración de punto fijo
Iteración de Newton
Conclusiones
2
/48
EL DOMINIO DE ESTABILIDAD NUMÉRICA
Basado en el libro de Cellier, F.E. and E. Kofman (2006), Continuous System Simulation, Springer-Verlag, New York
3
/48
Dominio de estabilidad analítica
4
Un sistema lineal invariante en el tiempo autónomo puede ser representado usando el modelo
La solución analítica es
La solución es analíticamente estable si
Región de estabilidad analítica en el plano λ
/485
•Dominio de estabilidad numérica para Euler directo
/48
Integración numérica usando Euler directo
6
Usando el algoritmo de Euler directo:
El sistema original continuo se ha convertido en un sistema de tiempo discreto “equivalente”
/48
Valores propios del sistema de tiempo discreto
7
El sistema continuo autonomo original es
El sistema de tiempo discreto autonomo “equivalente” es
Si λ es un valor propio de A, entonces
1 + λ·h es un valor propio de F
Demostrar!
/48
Sistema de tiempo discreto equivalente
8
En general, el sistema lineal de tiempo discreto “equivalente” puede expresarse usando
Región de estabilidad numerica en el plano λ·h
La region de estabilidad de un sistema de tiempo discreto es un circulo de radio unitario centrado en el origen
/48
Simulación con el algoritmo de Euler directo
9
1 h = 1
Simulacion para valores de a = - 0.1, -1, -2 y -3 con paso de integracion fijo h = 1
/48
Simulación con el algoritmo de Euler directo
10
Simulación con el algoritmo de Euler directo del sistema:
Con a = - 3.0 , a fin de obtener un resultado una exactitud del 0.15 %, el tamaño del paso debe ser h = 0.0001
¡10.000 pasos de integración para simular
10 segundos!
/48
El mayor paso de integración posible para Euler directo
11
Dado un sistema lineal de segundo orden con dos valores propios complejos λ1 y λ2
Región de estabilidad numerica en el plano λ·h
/48
El mayor paso de integración posible para Euler directo
12
Dado un sistema lineal con valores propios complejos λi :
Ejercicio
Demuestre que el tamaño del paso hmarg para el cual el
algoritmo directo de Euler dara marginalmente estable es:
Ayuda: usar el teorema de Tales
/48
Euler directo para sistemas marginalmente estables
13
La limitacion en el valor de h es particularmente importante cuando los autovalores se encuentran cerca del eje imaginario.
d se hace arbitrariamente pequeño
Cuando los autovalores estan sobre el eje imaginario no existe ningun paso de integracion que permita obtener una solucion
puramente oscilante
/4814
•Dominio de estabilidad numérica para Euler inverso
/48
Integración numérica usando Euler inverso
15
Usando el algoritmo de Euler inverso:
El sistema original continuo se ha convertido en un sistema de tiempo discreto “equivalente”
/48
Valores propios del sistema de tiempo discreto
16
El sistema continuo autonomo original es
El sistema de tiempo discreto autonomo “equivalente” es
¿Si λ es un valor propio de A, cual es el valor propio de F ?
/48
Región de estabilidad del Euler inverso
17
Dominio de estabilidad numérica del algoritmo de Euler inverso
Región de estabilidad numerica en el plano λ·h
El tamaño del paso de integración es dictado
sólo por requerimientos de exactitud,
No por estabilidad numérica
/48
Simulación con el algoritmo de Euler inverso
18
Simulación con el algoritmo de Euler inverso del sistema:
Este tipo de algoritmo es más apropiado que el Euler directo para resolver problemas con valores propios lejanos sobre el
real negativo en el plano λ
1 h = 1
/48
Euler inverso y los sistemas Stiff
19
Región de estabilidad numerica en el plano λ·h
Apropiado para sistemas stiff, es decir, sistemas con autovalores cuyas partes reales estan desparramadas a lo largo
del eje real negativo
/48
Simulación con el algoritmo de Euler inverso
20
Simulación con el algoritmo de Euler inverso del sistema:
Lección a aprender: Puede ser buena idea simular el sistema dos veces:
La simulación sugiere un sistema perfectamente
estable
Con a = +3.0
el sistema original es inestable
Una vez con algoritmo que exhiba un dominio de estabilidad comparable con el algoritmo de Euler directo,
Y otra con un algoritmo que se comporte como el algoritmo de Euler inverso
/48
IMPLEMENTACIÓN EN MATALB
21
•Cálculo del dominio de estabilidad numérica
/48
Definición de la Matrix del sistema
22
Considerando un sistema lineal de segundo orden con dos valores propios complejos λ1 y λ2 con la matriz
λ1 y λ2 estan localizados sobre el círculo unitariolambda = cos(α) + j*sen(α)
cos(α) – j*sen(α)
/48
Definición de la Matrix del sistema
23
Considerando un sistema lineal de segundo orden con dos valores propios complejos λ1 y λ2 con la matriz
Implementación en MATLAB
/48
Cálculo de la matriz F
24
Euler directo
Euler inverso
/48
Cálculo del máximo valor de h
25
Con h = hmax , los valores propios de F se encuentran en el
círculo unitario
Esta función trabaja para algoritmos con dominios de estabilidad similares al algoritmo Euler
/48
Cálculo del dominio de estabilidad numérica
26
Ejercicio:
Hacer un grafico del valor de hmax en función de α,
para el algoritmo Euler directo
en coordenadas polares
Utilice las funciones aa.m, ff.m, hh.m , y stabdom.m
/48
ITERACIÓN DE PUNTO FIJO
27
/48
Euler inverso para sistemas no lineales
28
En el caso de los sistemas lineales, en la simulación del algoritmo de Euler inverso la matriz F puede calcularse explícitamente
En el caso no lineal esto no puede hacerse
En el caso no lineal es necesario resolver un conjunto de ecuaciones algebraicas no lineales
implícitas
/48
Iteración de punto fijo
29
Una posible aproximación:
* Iniciar con una predicción
* Continuar con la iteración de varias correcciones
01 ( , )k k k kx x h f x t
predictor:
corrector:
11 1 1,i i
k k k kx x h f x t
Iteration i
/48
Iteración de punto fijo
30
Una posible aproximación:
/48
Iteración de punto fijo: caso lineal
31
Cuando se aplica la reiteración de punto fijo para un sistema lineal, se tiene
Después de un número infinito de interacciones
restando las dos ecuaciones
entonces
La misma matriz F que en el caso del algoritmo de Euler
inverso
/48
Iteración de punto fijo: caso lineal
32
Dominio de estabilidad numérica de la técnica predictor-corrector obtenido con el algoritmo que genera los dominios de estabilidad
/48
Iteración de punto fijo: caso lineal
33
La aproximación no funciona porque la serie infinita
Dentro del círculo unitario, el dominio de estabilidad numérica del método predictor-corrector es el mismo que para algoritmo
de Euler inverso
Fuera del círculo unitario, el método es inestable
sólo converge si todos los valores propios de A·h se encuentran dentro del círculo unitario
Si no es el caso, la resta de las ecuaciones es inválida
/48
ITERACIÓN DE NEWTON
34
/48
Iteración de Newton
35
La iteración de Newton puede ser usada para determinar el cruce por cero de una función
Cruce por cero
/48
Iteración de Newton: caso escalar
36
El algoritmo de Euler inverso aplicado a una ecuación diferencial escalar es
Por lo tanto:
Aplicando la iteración de Newton
donde k es el numero del paso de integracion y l es el numero de veces que la iteracion de Newton fue aplicada en dicho paso
/48
Iteración de Newton: extensión matricial
37
La extensión matricial del algoritmo de iteración de Newton es
Dada en términos de la matriz Hessiana de la iteración de Newton
/48
Iteración de Newton: caso vectorial
38
Aplicando la iteración de Euler al vector del modelo en espacio de estado y al algoritmo de Euler inverso, se tiene
Dada en términos de la matriz Jacobiana del sistema dinámico
/48
Iteración de Newton: caso vectorial
39
Matriz Hessiana
Jacobiana
Aplicando la iteración de Euler al vector del modelo en espacio de estado y al algoritmo de Euler inverso
inversión
En cada iteración de Euler, en general, es necesario calcular la matriz Jacobiana, así como invertir la matriz Hessiana
/48
Iteración de Newton: caso vectorial
40
Si el sistema es lineal
La iteración de Newton no cambia el dominio de estabilidad numérica de un solver.
Esto es cierto no solamente para el algoritmo de Euler inverso, si no tambien para todos los solvers numéricos
Entonces:
/48
CONCLUSIONES
41
/48
Conclusiones
42
En el análisis de un solver numérico, la estabilidad numérica del algoritmo debe ser siempre tomada en consideración
La estabilidad numérica de la mayoría de los solvers puede ser representada por un dominio de estabilidad numérica, dibujado en el plano complejo λ·h
/48
Conclusiones
43
En el análisis de un solver numérico, la estabilidad numérica del algoritmo debe ser siempre tomada en consideración
La estabilidad numérica de la mayoría de los solvers puede ser representada por un dominio de estabilidad numérica, dibujado en el plano complejo λ·h
La estabilidad numérica de los solvers es usualmente analizada solamente para sistemas lineales autónomos invariantes en el tiempo
/48
Conclusiones
44
En el análisis de un solver numérico, la estabilidad numérica del algoritmo debe ser siempre tomada en consideración
/48
Conclusiones
45
En el análisis de un solver numérico, es también importante considerar la exactitud de la aproximación del algoritmo
La exactitud numérica de un solver p está sujeta a diferentes tipos de error, tales como el error por truncado, el error por redondeo, y el error por acumulación
/48
Conclusiones
46
En el análisis de un solver numérico, es también importante considerar la exactitud de la aproximación del algoritmo
La exactitud numérica de un solver p está sujeta a diferentes tipos de error, tales como el error por truncado, el error por redondeo, y el error por acumulación
El más importante de estos errores es el error por truncado, que está caracterizado por el orden de la exactitud de la aproximación del solver
/48
Fuentes
Cellier, F.E. and E. Kofman (2006), Continuous System Simulation, Springer-Verlag, New York
47
/48
FIN
48