PDF de programación - Un caso de modelado y simulación

Un caso de modelado y simulación

F. Javier Gil Chica

abril, 2009

Resumen

Este artículo se ofrece como introducción a la asignatura de Mode-
lado y simulación de sistemas dinámicos. Un sistema dinámico es aquél
que se describe mediante un cierto número de variables cuyos valores
son función del tiempo. Las teorías que se ocupan de diversas clases
de sistemas dinámicos (así, la Mecánica para sistemas mecánicos, la
Electricidad para sistemas eléctricos, etc.) no nos proporcionan direc-
tamente estas funciones, sino ecuaciones que expresan la velocidad a
la que cambian las variables que describen el sistema. Formalmente,
nos proporcionan sistemas de ecuaciones diferenciales que después han
de ser integradas para obtener propiamente las funciones que expre-
san los valores de las variables en función del tiempo. A veces, esta
integración puede hacerse de forma analítica. A veces, ha de hacerse
de forma numérica.

1. Nuestro sistema

Elegimos como sistema ilustrativo un proyectil de masa m. Desea-
mos conocer su posición y velocidad en todo instante, y quizás algún
dato derivado, como su máximo alcance o máxima altura alcanzada.
El proyectil puede modelarse de diversas formas (como una masa iner-
te lanzada con una velocidad inicial o como proyectil autopropulsado)
y ponerse en diversos escenarios: movimiento ideal sin rozamiento y
sometido a la aceleración constante de la gravedad, con rozamiento
dependiente de la velocidad, con rozamiento que depende de la velo-
cidad y de la densidad del aire, con una densidad del aire que varía de
una forma u otra según la altura, lanzado desde una plataforma fija o
desde una plataforma en movimiento, sometido o no a los efectos de
la rotación del sistema de referencia, etc.

1

Dependiendo del modelo de proyectil elegido y del escenario las
ecuaciones que describen su movimiento tendrán una u otra forma. Si
las ecuaciones son sencillas, será fácil integrarlas. En caso contrario,
será preciso integrarlas numéricamente. Como resultado de la inte-
gración numérica, obtenemos listas de números que hemos de com-
binar para obtener gráficas. La integración numérica se efectúa me-
diante un programa. Aquí usaremos el lenguaje C y representaremos
las gráficas mediante gnuplot. Hay otras opciones: Matlab, Maple,
Mathematica...

El procedimiento por tanto será el siguiente:

1. Identificar exactamente qué sistema queremos modelar y simular.

2. Con ayuda de la ciencia que se ocupe de los sistemas del tipo
al que pertenezca aquél en que estamos interesados, obtendre-
mos un modelo para su comportamiento. Si esto no es suficiente,
tendremos que hacer hipótesis razonables.

3. En general, cuando decimos modelo nos estamos refiriendo a un

sistema de ecuaciones diferenciales.

4. Si el modelo puede resolverse analíticamente, lo resolveremos así.
Si no, lo integraremos numéricamente y representaremos gráfi-
camente el resultado.

Hay otros muchos aspectos que pueden tratarse. Por ejemplo, po-
demos preguntarnos por la estabilidad del sistema, es decir, por su
comportamiento cuando t → ∞. Es evidente que no podemos integrar
numéricamente hasta t → ∞ y lo normal es que tampoco sepamos
integrar analíticamente el modelo, por lo que se requerirán técnicas
especiales que se tratan durante el curso. También es posible que no
nos baste conocer cómo se va a comportar el sistema, sino que desee-
mos manipularlo de alguna manera para que se comporte de alguna
forma útil para nosotros. Existen igualmente técnicas potentes para
conseguir esto que se tratarán también durante el curso.

2. Un método de integración

El sistema dinámico más sencillo que podamos concebir es aquél que
se describe mediante una única variable x. Con ayuda de hipótesis o
mediante alguna ciencia particular (Mecánica, Electricidad, Química,
etc.), obtendremos un modelo dinámico que se expresará como una
ecuación diferencial de la forma

2

˙x = f (x, t)

(1)

˙x indica derivada de x respecto al tiempo. La función f (x, t) puede
ser más o menos sencilla, lineal o no lineal, dependiente de x y t o de
sólo una de ellas, etc. El problema que se plantea es obtener la función
x(t). Existen multitud de algoritmos para la resolución de este tipo de
ecuaciones. Nosotros presentaremos uno de ellos, el método de Runge-
Kutta de cuarto orden. Es un método popular, sencillo de programar,
relativamente rápido y en la mayor parte de las ocasiones de suficiente
precisión. No vamos a discutir la forma en que se obtiene dicho método,
sino simplemente describir el algoritmo. Antes de eso, conviene aclarar
la naturaleza de todos estos métodos. Una integración numérica no
nos proporciona la función x(t), sino una serie de valores de x para
instantes t1, t2, t3...tn. Aunque existen métodos de paso variable donde
los intervalos ti+1 − ti dependen de i, en la mayoría de las ocasiones,
para todos los i, la diferencia ti+1 − ti es una constante h, de manera
que si partimos desde un instante inicial t0, ti = t0 + ih. Todos estos
métodos requieren el conocimiento de un estado inicial a partir del
cual calcular estados sucesivos. Particularizando en nuestro sistema,
es obvio que la posición y velocidad del proyectil en un instante dado
depende de la posición y velocidad con que fue lanzado.

El método de Runge-Kutta parte de un valor xk (indicamos así x(tk))
y calcula el valor xk+1. Si fijamos h (por ejemplo, h = 0,1 cuando de-
seamos conocer x de décima en décima de segundo) el algoritmo sigue
los siguientes pasos:

1. Calcular k1 = hf (xk, tk)

2. Calcular k2 = hf (xk + k1/2, tk + h/2)

3. Calcular k3 = hf (xk + k2/2, tk + h/2)

4. Calcular k4 = hf (xk + k3, tk + h)
5. Calcular xk+1 = xk + 1

6 [k1 + 2k2 + 2k3 + k4]

En general, los sistemas dinámicos vienen descritos no por una sola
variable, sino por varias. Por consiguiente, necesitamos generalizar el
algoritmo para el caso en que queremos integrar no una ecuación sino
un conjunto de ellas. Daremos el algoritmo explícito para un sistema
de dos ecuaciones. La generalización para un número cualquiera es
inmediata y se deja como ejercicio al lector. Sea el sistema:

3

˙x = f (x, y, t)

˙y = g(x, y, t)

(2)

Dados unos valores xk e yk, el algoritmo sigue los siguientes pasos:

1. Calcular k1 = hf (xk, yk, tk)

2. Calcular n1 = hg(xk, yk, tk)

3. Calcular k2 = hf (xk + k1/2, yk + n1/2, tk + h/2)

4. Calcular n2 = hg(xk + k1/2, yk + n1/2, tk + h/2)

5. Calcular k3 = hf (xk + k2/2, yk + n2/2, tk + h/2)

6. Calcular n3 = hg(xk + k2/2, yk + n2/2, tk + h/2)

7. Calcular k4 = hf (xk + k3, yk + n3, tk + h)

8. Calcular n4 = hg(xk + k3, yk + n3, tk + h)
9. Calcular xk+1 = xk + 1
10. Calcular yk+1 = yk + 1

6 [k1 + 2k2 + 2k3 + k4]
6 [n1 + 2n2 + 2n3 + n4]

3. Variables de estado

Sucede con frecuencia que nos disponemos de modelos que se expre-
sen como sistemas de ecuaciones diferenciales de primer orden, sino
de orden dos o superiores. Es el caso de los problemas de la mecánica.
Por ejemplo, de la ecuación de Newton:

¯F = m¯a

Escribiéndola en sus componentes:

¨x =

¨y =

¨z =

1
m
1
m
1
m

Fx

Fy

Fz

4

(3)

(4)

Este es un sistema de ecuaciones diferenciales de segundo orden,
ya que aparecen las segundas derivadas de las variables. Un sistema
de este tipo puede siempre reducirse a otro sistema de ecuaciones de
primer orden introduciendo variables de estado. Sea por ejemplo un
modelo de la forma ¨x = f (x, t). Si introducimos variables x1 y x2 de
forma que

es claro entonces que

x1 = x
x2 = ˙x

˙x1 = x2
˙x2 = f (x1, t)

(5)

(6)

Otro ejemplo: para un modelo de la forma x(3) = f (x, t), donde x(3)
representa la tercera derivada respecto al tiempo de x, introducimos
tres variables de estado: x1 = x, x2 = ˙x y x3 = ¨x, con lo cual:

˙x1 = x2
˙x2 = x3
˙x3 = f (x1, t)

(7)

Hay un caso particular interesante, y es cuando el modelo consiste

en una única ecuación diferencial ordinaria de grado n de la forma

dnx
dtn + an−1

dn−1x
dtn−1 + an−2

dn−2x
dtn−2 + · · · + a1

dx
dt

+ a0x = 0

(8)

Si introducimos variables de estado

x1 = x
x2 = ˙x
x3 = ¨x
... = ...

xn =

dn−1x
dtn−1

5

(9)

es claro que

˙x1 = x2
˙x2 = x3
˙x3 = x4
... = ...
˙xn = −a0x1 − a1x2 − a2x3 − · · · − an−1xn

(10)

4. Primer modelo

Nuestro primer modelo para el proyectil es muy simple: una masa m
que se lanza con una velocidad inicial desde el origen de coordenadas,
sin rozamiento con el aire, en un campo de gravedad constante. Puesto
que la gravedad actúa sólo en la dirección del eje z, el movimiento
resultante está contenido en el plano (x, z). El proyectil es lanzado
en la dirección positiva del eje x y el eje z está dirigido verticalmente
hacia arriba. La ecuación de Newton, resuelta en sus componentes, nos
dice que la aceleración según el eje x es nula, y que la aceleración según
el eje z es −g, siendo g el módulo de la aceleración de la gravedad:

¨x = 0

¨z = −g

(11)

Necesitamos introducir las varibles de estado: x1 = x, x2 = z, x3 =

˙x y x4 = ˙z, con lo cual

˙x1 = x3
˙x2 = x4
˙x3 = 0
˙x4 = −g

Si llamamos ¯x al vector de estado de componentes

¯x =

x1
x2
x3
x4







6

(12)

(13)

el sistema anterior puede escribir en la forma matricial compacta

siguiente:

˙¯x =




0 0 1 0
0 0 0 1
0 0 0 0
0 0 0 0




¯x +

0
0
0
−g







(14)

Este es un caso particular de un tipo de sistemas de la forma general

˙¯x = A¯x + B¯u

(15)

donde A y B son ciertas matrices. Existe una solución analítica pa-
ra los sistemas de este tipo y de ella se tratará durante el curso. Este
sistema se integra de forma trivial, de manera que podemos obtener
expresiones analíticas para x(t) e y(t), y de ahí obtener algunos resul-
tados de interés, como el máximo alcance, altura máxima, ecuación
de la envolvente, tiempo de vuelo, etc. Estos resultados se encuentran
en los textos elementales, así que no los discutiremos aquí. En lugar
de eso, implementaremos un programa que integre numéricamente el
modelo