Un caso de modelado y simulaci on

F. Javier Gil Chica abril, 2009

Resumen Este art culo se ofrece como introducci on a la asignatura de Modelado y simulaci on de sistemas din amicos. Un sistema din amico es aqu el que se describe mediante un cierto n umero de variables cuyos valores son funci on del tiempo. Las teor as que se ocupan de diversas clases de sistemas din amicos (as , la Mec anica para sistemas mec anicos, la Electricidad para sistemas el ectricos, etc.) no nos proporcionan directamente 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 es han de ser integradas para obtener propiamente las funciones que expresan los valores de las variables en funci on del tiempo. A veces, esta integraci on puede hacerse de forma anal tica. A veces, ha de hacerse de forma num erica.

1.

Nuestro sistema

Elegimos como sistema ilustrativo un proyectil de masa m. Deseamos conocer su posici on y velocidad en todo instante, y quiz as alg un dato derivado, como su m aximo alcance o m axima altura alcanzada. El proyectil puede modelarse de diversas formas (como una masa inerte lanzada con una velocidad inicial o como proyectil autopropulsado) y ponerse en diversos escenarios: movimiento ideal sin rozamiento y sometido a la aceleraci on constante de la gravedad, con rozamiento dependiente de la velocidad, con rozamiento que depende de la velocidad y de la densidad del aire, con una densidad del aire que var a de una forma u otra seg un la altura, lanzado desde una plataforma ja o desde una plataforma en movimiento, sometido o no a los efectos de la rotaci on del sistema de referencia, etc.

Dependiendo del modelo de proyectil elegido y del escenario las ecuaciones que describen su movimiento tendr an una u otra forma. Si las ecuaciones son sencillas, ser a f acil integrarlas. En caso contrario, ser a preciso integrarlas num ericamente. Como resultado de la integraci on num erica, obtenemos listas de n umeros que hemos de combinar para obtener gr acas. La integraci on num erica se efect ua mediante un programa. Aqu usaremos el lenguaje C y representaremos las gr acas mediante gnuplot. Hay otras opciones: Matlab, Maple, Mathematica... El procedimiento por tanto ser a el siguiente: 1. Identicar exactamente qu e sistema queremos modelar y simular. 2. Con ayuda de la ciencia que se ocupe de los sistemas del tipo al que pertenezca aqu el en que estamos interesados, obtendremos un modelo para su comportamiento. Si esto no es suciente, tendremos que hacer hip otesis razonables. 3. En general, cuando decimos modelo nos estamos reriendo a un sistema de ecuaciones diferenciales. 4. Si el modelo puede resolverse anal ticamente, lo resolveremos as . Si no, lo integraremos num ericamente y representaremos gr acamente el resultado. Hay otros muchos aspectos que pueden tratarse. Por ejemplo, podemos preguntarnos por la estabilidad del sistema, es decir, por su comportamiento cuando t . Es evidente que no podemos integrar num ericamente hasta t y lo normal es que tampoco sepamos integrar anal ticamente el modelo, por lo que se requerir an t ecnicas especiales que se tratan durante el curso. Tambi en es posible que no nos baste conocer c omo se va a comportar el sistema, sino que deseemos manipularlo de alguna manera para que se comporte de alguna forma u til para nosotros. Existen igualmente t ecnicas potentes para conseguir esto que se tratar an tambi en durante el curso.

2.

Un m etodo de integraci on

El sistema din amico m as sencillo que podamos concebir es aqu el que se describe mediante una u nica variable x. Con ayuda de hip otesis o mediante alguna ciencia particular (Mec anica, Electricidad, Qu mica, etc.), obtendremos un modelo din amico que se expresar a como una ecuaci on diferencial de la forma

x = f (x, t)

(1)

x indica derivada de x respecto al tiempo. La funci on f (x, t) puede ser m as o menos sencilla, lineal o no lineal, dependiente de x y t o de s olo una de ellas, etc. El problema que se plantea es obtener la funci on x(t). Existen multitud de algoritmos para la resoluci on de este tipo de ecuaciones. Nosotros presentaremos uno de ellos, el m etodo de RungeKutta de cuarto orden. Es un m etodo popular, sencillo de programar, relativamente r apido y en la mayor parte de las ocasiones de suciente precisi on. No vamos a discutir la forma en que se obtiene dicho m etodo, sino simplemente describir el algoritmo. Antes de eso, conviene aclarar la naturaleza de todos estos m etodos. Una integraci on num erica no nos proporciona la funci on x(t), sino una serie de valores de x para instantes t1 , t2 , t3 ...tn . Aunque existen m etodos 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 etodos requieren el conocimiento de un estado inicial a partir del cual calcular estados sucesivos. Particularizando en nuestro sistema, es obvio que la posici on y velocidad del proyectil en un instante dado depende de la posici on y velocidad con que fue lanzado. El m etodo de Runge-Kutta parte de un valor xk (indicamos as x(tk )) y calcula el valor xk+1 . Si jamos h (por ejemplo, h = 0,1 cuando deseamos conocer x de d ecima en d ecima 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)
1 [k1 + 2k2 + 2k3 + k4 ] 5. Calcular xk+1 = xk + 6

En general, los sistemas din amicos 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 on sino un conjunto de ellas. Daremos el algoritmo expl cito para un sistema de dos ecuaciones. La generalizaci on para un n umero cualquiera es inmediata y se deja como ejercicio al lector. Sea el sistema:

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)
1 [k1 + 2k2 + 2k3 + k4 ] 9. Calcular xk+1 = xk + 6

10. Calcular yk+1 = yk + 1 6 [n1 + 2n2 + 2n3 + n4 ]

3.

Variables de estado

Sucede con frecuencia que nos disponemos de modelos que se expresen como sistemas de ecuaciones diferenciales de primer orden, sino de orden dos o superiores. Es el caso de los problemas de la mec anica. Por ejemplo, de la ecuaci on de Newton: = ma F Escribi endola en sus componentes: 1 Fx m 1 y = Fy m 1 z = Fz m (3)

x =

(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 x1 = x x2 = x es claro entonces que x 1 = x2 x 2 = f (x1 , t) (6) (5)

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: x 1 = x2 x 2 = x3 x 3 = f (x1 , t) (7)

Hay un caso particular interesante, y es cuando el modelo consiste en una u nica ecuaci on diferencial ordinaria de grado n de la forma dn1 x dn2 x dx dn x + a + a + + a1 + a0 x = 0 n 1 n 2 n n 1 n 2 dt dt dt dt Si introducimos variables de estado x1 = x x2 = x x3 = x ... = ... dn1 x xn = dtn1 (9)

(8)

es claro que x 1 = x2 x 2 = x3 x 3 = x4 ... = ... x n = a0 x1 a1 x2 a2 x3 an1 xn (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 ua s olo en la direcci on del eje z , el movimiento resultante est a contenido en el plano (x, z ). El proyectil es lanzado en la direcci on positiva del eje x y el eje z est a dirigido verticalmente hacia arriba. La ecuaci on de Newton, resuelta en sus componentes, nos dice que la aceleraci on seg un el eje x es nula, y que la aceleraci on seg un el eje z es g, siendo g el m odulo de la aceleraci on 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 x 1 = x3 x 2 = x4 x 3 = 0 x 4 = g Si llamamos x al vector de estado de componentes

(12)

x =

x1 x2 x3 x4

(13)

el sistema anterior puede escribir en la forma matricial compacta siguiente:

= x

0 0 0 0

0 0 0 0

1 0 0 0

0 1 0 0

+ x

0 0 0 g

(14)

Este es un caso particular de un tipo de sistemas de la forma general = Ax x + Bu (15)

donde A y B son ciertas matrices. Existe una soluci on anal tica para los sistemas de este tipo y de ella se tratar a 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 resultados de inter es, como el m aximo alcance, altura m axima, ecuaci on 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 ericamente el modelo. Seguiremos los pasos desde la resoluci on num erica hasta la obtenci on de una gr aca. Despu es, usaremos esa gr aca para comparar este primer modelo con modelos m as elaborados. Puesto que con frecuencia la visualizaci on de los datos es relevante, nos detendremos en los detalles del proceso que lleva desde la resoluci on num erica hasta la generaci on del documento nal que incluye una gr aca. Por supuesto, existen muchas formas de hacer esto, y seguramente m as c omodas que la que adoptamos aqu . El siguiente es el listado del programa que implementa el algoritmo de Runge-Kutta para resolver el sistema de ecuaciones (12) 1 .

Por conveniencia, hemos suprimido la indentaci on original, para que el c odigo fuente pueda adaptarse al ancho de este documento. No ser a un grave inconveniente, puesto que el programa es corto

#include <stdio.h> #include <math.h> #define pi 3.14159265 #define G 9.81 #define V0 100 /* aceleracion gravedad */ /* velocidad inicial */

double k[4],l[4],m[4],n[4]; double h=0.01; double t=0; double x0[4], x1[4]; FILE *listado;

double f0(double x0, double x1, double x2, double x3, double t) { return(x2); } double f1(double x0, double x1, double x2, double x3, double t) { return(x3); } double f2(double x0, double x1, double x2, double x3, double t) { return(0); } double f3(double x0, double x1, double x2, double x3, double t) { return(-G); } main() { if ((listado=fopen("datos.dat","w"))==NULL) return; /* establecer condiciones iniciales */ x0[0]=0; x0[1]=0; x0[2]=V0*cos(pi/4); x0[3]=V0*sin(pi/4);

/* mientras la coordenada z sea >=0 */ while (x0[1]>=0){ /* salvar a archivo */ fprintf(listado,"%f %f %f %f %f\n",t,x0[0],x0[1],x0[2],x0[3]); /* paso 1 */ k[0]=h*f0(x0[0],x0[1],x0[2],x0[3],t); l[0]=h*f1(x0[0],x0[1],x0[2],x0[3],t); m[0]=h*f2(x0[0],x0[1],x0[2],x0[3],t); n[0]=h*f3(x0[0],x0[1],x0[2],x0[3],t); /* paso 2 */ k[1]=h*f0(x0[0]+k[0]/2,x0[1]+l[0]/2,x0[2]+m[0]/2,x0[3]+n[0]/2,t+h/2); l[1]=h*f1(x0[0]+k[0]/2,x0[1]+l[0]/2,x0[2]+m[0]/2,x0[3]+n[0]/2,t+h/2); m[1]=h*f2(x0[0]+k[0]/2,x0[1]+l[0]/2,x0[2]+m[0]/2,x0[3]+n[0]/2,t+h/2); n[1]=h*f3(x0[0]+k[0]/2,x0[1]+l[0]/2,x0[2]+m[0]/2,x0[3]+n[0]/2,t+h/2); /* paso 3 */ k[2]=h*f0(x0[0]+k[1]/2,x0[1]+l[1]/2,x0[2]+m[1]/2,x0[3]+n[1]/2,t+h/2); l[2]=h*f1(x0[0]+k[1]/2,x0[1]+l[1]/2,x0[2]+m[1]/2,x0[3]+n[1]/2,t+h/2); m[2]=h*f2(x0[0]+k[1]/2,x0[1]+l[1]/2,x0[2]+m[1]/2,x0[3]+n[1]/2,t+h/2); n[2]=h*f3(x0[0]+k[1]/2,x0[1]+l[1]/2,x0[2]+m[1]/2,x0[3]+n[1]/2,t+h/2); /* paso 4 */ k[3]=h*f0(x0[0]+k[2],x0[1]+l[2],x0[2]+m[2],x0[3]+n[2],t+h); l[3]=h*f1(x0[0]+k[2],x0[1]+l[2],x0[2]+m[2],x0[3]+n[2],t+h); m[3]=h*f2(x0[0]+k[2],x0[1]+l[2],x0[2]+m[2],x0[3]+n[2],t+h); n[3]=h*f3(x0[0]+k[2],x0[1]+l[2],x0[2]+m[2],x0[3]+n[2],t+h); /* nuevos valores */ x1[0]=x0[0]+(1.0/6.0)*(k[0]+2*(k[1]+k[2])+k[3]); x1[1]=x0[1]+(1.0/6.0)*(l[0]+2*(l[1]+l[2])+l[3]); x1[2]=x0[2]+(1.0/6.0)*(m[0]+2*(m[1]+m[2])+m[3]); x1[3]=x0[3]+(1.0/6.0)*(n[0]+2*(n[1]+n[2])+n[3]); /* actualizar */ x0[0]=x1[0]; x0[1]=x1[1]; x0[2]=x1[2]; x0[3]=x1[3]; t=t+h; } /* end while */ fclose(listado); return; }

Hemos tomado un incremento h = 0,01. El proyectil es lanzado desde el origen de coordenadas con una velocidad de 100 m/s y formando un angulo de /4 con el eje horizontal. El programa calcula mientras el proyectil se encuentre sobre el suelo, es decir, mientras sea z >= 0. Los datos se guardan en el archivo datos.dat, que es un archivo de texto que contiene cinco columnas, a saber: el tiempo, la coordenada x, la coordenada z , la velocidad seg un el eje x y la velocidad seg un el eje z . Hemos usado el programa gnuplot para representar gr acamente la trayectoria del proyectil, es decir, para representar la gr aca z (x). La gr aca puede representarse en pantalla (la opci on por defecto) pero puede tambi en guardarse en otros formatos. Hemos elegido PostScript encapsulado porque podemos incorporar esos gr acos a nuestro documento .tex mediante el paquete dvips. Una vez dentro de gnuplot, escribimos > > > > > > set terminal fig set output "trayectoria-0.fig" set xlabel "x" set ylabel "z(x)" set grid plot "datos.dat" using 2:3 with dots notitle

La primera l nea indica al programa que la salida no ser a a pantalla, sino en el formato propio del programa xfig. La segunda l nea establece el nombre del archivo gr aco. Las l neas tercera y cuarta establecen las etiquetas para los ejes horizontal y vertical. La quinta l nea fuerza a que se dibuje la cuadr cula, que por defecto se omite. Una vez hecho esto, la u ltima l nea es la que se encarga efectivamente de trazar la gr aca, mediante la orden plot, que toma los par ametros siguientes: datos.dat es el nombre del archivo donde se encuentran los datos; using 2:3 indica a plot que, de todas las columnas que componen el archivo, use la segunda y la tercera; la segunda contiene los valores de x y la tercera los valores de z ; with dots indica a plot que no trace una l nea continua, uniendo puntos, sino que cada pareja (x, z ) sea representada por un punto; notitle indica a plot que omita en la gr aca una etiqueta que aparece por defecto con el nombre del archivo de donde provienen los datos. Una vez generado el archivo gr aco, podemos editarlo y modicarlo si fuese preciso desde xfig y a continuaci on lo exportamos como archivo .eps y procedemos a incluirlo en nuestro documento fuente .tex. Para ello, introducimos las siguientes l neas:

10

\begin{figure}{h} \begin{center} \includegraphics[scale=1.0]{trayectoria-0.eps}\\ {\bf Figura 1} \end{center} \end{figure} El resultado es el que aparece:
300

250

200

z(x)

150

100

50

200

400

600

800

1000

1200

Figura 1

5.

Segundo modelo

Modelaremos el movimiento del proyectil suponiendo que existe un rozamiento proporcional al cuadrado del m odulo de la velocidad y de sentido contrario a la misma. Por consiguiente, la ecuaci on de Newton incluye tanto la fuerza de la gravedad como la resistencia del medio: ma = mg kv v (16)

Simplicando la masa, introduciendo la constante k = k/m y escribiendo expresiones expl citas para las componentes x y z tenemos que x = k v x z = k v z (17)

11

donde |v | = x 2 + z 2 (18)

Introduciendo variables de estado x1 = x, x2 = z , x3 = x y x4 = z : x 1 = x3 x 2 = x4 x 3 = k vx3 x 4 = g k vx4 (19)

Obs ervese que, respecto al programa presentado anteriormente, s olo hemos de modicar el cuerpo de las funciones f0(), f1(), f2() y f3(). El resto queda igual. En la Figura 2 puede apreciarse el efecto del rozamiento sobre el alcance y la m axima altura, para los valores k = 0 (caso sin rozamiento), k = 0,0001 y k = 0,001:
300

250

k=0 k=104

200

z(x)

150

k=103

100

50

200

400

600

800

1000

1200

Figura 2

6.

Tercer modelo

Hemos considerado en el modelo anterior que el movimiento del proyectil sufre una fuerza de rozamiento que depende del cuadrado de la velocidad a la que se mueve. Tambi en depende de la densidad del medio en el que se mueve. Si esta densidad es constante, podemos

12

considerarla integrada en la constante k , de forma que k = c, siendo c otra constante. Pero sabemos que la densidad del aire var a con la altura, de forma que k = c(z ). Hemos de buscar una forma para la funci on (z ). Podr amos acudir a la hidrost atica pero, como no deseamos entrar en ese punto, haremos alguna hip otesis razonable. Sabemos que la densidad es m axima en la supercie y que decrece progresivamente (no hay borde en la atm osfera) hasta que desaparece pr acticamente a unos 100 km (el l mite es algo arbitrario). Una funci on que se adapta a estas caracter sticas es una exponencial: (z ) = (0) exp(z ) (20)

Queda pendiente el ajuste de las constantes. Para ajustar , usaremos el hecho de que la densidad del aire se reduce a la mitad aproximadamente a unos 5000 m. de altura: (5000) = de donde = 1 1 ln( ) 1,386 104 5000 2 (22) 1 (0) = (0) exp(5000) 2 (21)

Para ajustar c, podemos hacerlo conjuntamente con (0), de forma que k = c(0) y coincida num ericamente con el valor de k usado en el modelo anterior. De esta forma, la simulaci on de este modelo requiere unicamente la sustituci on de k por k exp(z ). Por otro lado, el efecto ser a tanto m as evidente cuanta mayor sea la altura alcanzada por el proyectil. Para conseguir mayor altura, elevamos el a ngulo de tiro desde /4 a 3/8 e incrementamos la velocidad desde los 100 m/s hasta los 1000 m/s. Repetimos la simulaci on del segundo modelo con los nuevos valores de angulo y velocidad y comparamos con el modelo actual. El resultado de la simulaci on se presenta en la Figura 3

7.

Cuarto modelo

Nos planteamos ahora el movimiento de un proyectil autopropulsado. Un proyectil de tal clase ha de llevar consigo el combustible. Ese combustible se consume a un ritmo determinado y el empuje resultante es proporcional a ese ritmo. A medida que se consume el combustible, disminuye la masa del proyectil, de forma que ya no puede considerarse esta constante. El modelo por tanto ha de partir de

13

2500

densidad decreciente
2000

1500

z(x)

densidad constante
1000

500

200

400

600

800

1000

1200

1400

1600

Figura 3

alguna formulaci on m as general que la forma b asica de la ecuaci on de = ma Newton F . No es preciso ir muy lejos: = dp F dt (23)

con p = m(t) v . Si el ritmo al que se quema combustible es constante: m(t) = m(0) t tenemos = v p + m(t)v (25) (24)

Como fuerzas, consideraremos la fuerza de la gravedad, la fuerza de rozamiento proporcional al cuadrado de la velocidad y la fuerza de autopropulsi on, tangente a la trayectoria en todo punto: v =g F k vv + v Y este es el modelo: g k vv + v = v + m(t)v v (27) (26)

14

Introduciremos algunos par ametros descriptivos. En primer lugar, si la masa inicial es m(0), llamaremos f a la fracci on de esa masa que es combustible, y al tiempo durante el cual ese combustible va a proporcionar empuje. En cuanto a ese empuje, lo consideraremos constante. De aqu : m( ) = (1 f )m(0) = m(0) de donde = f m(0) (29) (28)

Por otra parte, el m odulo del empuje ha de ser superior al peso del 2 proyectil , luego > m(0)g de donde > g f (31) (30)

Introduciendo el factor p > 1 escribimos: =p g f (32)

Con estas deniciones, escribimos el modelo en la forma: g kv v + pm(0)g f f v = m(0) v + m(0)(1 t)v v (33)

Que se desdobla en las componentes: x v z m(0)gt kv z + pm(0)g v kv x + pm(0)g f = m(0)x + m(0)t x f = m(0)z + m(0)t z

(34)

donde por aligerar la escritura hemos hecho t = 1

2

f t

(35)

En rigor, no es as , pero prescindimos de los detalles

15

Introduciendo variables de estado tenemos: x 1 = x3 x 2 = x4 x 3 = x 4 kv pg f + + x3 m(0)t vt t pg f kv + + x4 = g + m(0)t vt t

(36)

Hemos realizado algunas simulaciones de este sistema, con valores m = 100kg., p = 3, f = 0,6, = 200s y v (0) = 10m/s. De nuevo, el programa que presentamos en el primer modelo s olo precisa modicaciones en las deniciones de las funciones f0(), f1(), f2() y f3(). Hay algunos aspectos cualitativos nuevos. En primer lugar, es preciso comprobar en cada paso de integraci on si queda a un combustible o no. El combustible se agota cuando t = . En segundo lugar, el proyectil no es disparado, sino que abandona una rampa de lanzamiento por sus propios medios. La velocidad a la salida de la rampa es comparativamente baja, y por tanto comparativamente a un proyectil no autopropulsado del mismo alcance el angulo de salida ha de ser m as elevado. La Figura 4 presenta la trayectoria para dos angulos distintos: 45o y 50o .
60

50

40

=50

z(x)

30

=45
20

10

100

200

300

400

500

600

700

800

900

Figura 4

Vemos que el alcance es sensible respecto al angulo de salida. De hecho, esa diferencia de 5o pr acticamente duplica el alcance. M as to-

16

dav a, con un angulo de 80o el alcance sube hasta los 210 km. y la altura m axima que se alcanza es de aproximadamente 21000 m. En tercer lugar, para la altura alcanzada es ya muy apreciable la disminuci on de la densidad del aire, y por tanto de la resistencia. Este factor no lo hemos tenido en cuenta en esta simulaci on. En cuarto lugar, para un alcance de m as de 200 km. ya es apreciable el efecto de la curvatura de la Tierra, y de su rotaci on. En quinto lugar, es preciso examinar el tiempo de vuelo. Para los ngulos de salida de 45o y 50o estos tiempos son aproximadamente de a 6 y 8 s., y eso signica que s olo una peque na fracci on del combustible se ha gastado, lo cual no tiene sentido. Para el angulo m as elevado de 80o el tiempo de vuelo es de 160 s., quedando a un una fracci on no despreciable de combustible por consumir. Efectuamos entonces otra simulaci on con un angulo de salida de 85o y comprobamos c omo el alcance y la altura m axima se incrementan, tal y como se muestra en la Figura 5:
120000

100000

80000

=85

z(x)

60000

40000

20000

50000

100000

150000

200000

250000

300000

350000

Figura 5

Es tambi en interesante examinar la curva x(t), que se muestra en la Figura 6:

17

350000 300000 250000 200000 150000 100000 50000 0

x(t)

50

100

150

200

250

300

350

400

450

Figura 6

Como se ve, el tiempo de vuelo supera los 400 s. y la curva es c oncava hasta t = 200s., indicando la aceleraci on del motor, pasando a convexa a partir de ese punto, indicando la deceleraci on por rozamiento. Este efecto se aprecia con mucha mayor claridad si se examina la gr aca para la velocidad en funci on del tiempo, representada en la Figura 7:
2000 1800 1600 1400 1200

x(t)

1000 800 600 400 200 0 0 50 100 150 200 250 300 350 400 450

Figura 7

Tambi en es interesante examinar la gr aca para la velocidad vertical, Figura 8:

18

1000 800 600 400

z(t)

200 0 200 400 600 800 0 50 100 150 200 250 300 350 400 450

Figura 8

8.

Quinto modelo

Las ecuaciones de Newton son v alidas para sistemas de referencia inerciales. Pero cuando un proyectil es lanzado desde tierra, lo hace desde una plataforma giratoria, sometida por tanto a aceleraci on, es decir, desde un sistema no inercial en el que no son v alidas las ecuaciones de Newton. Sea (x, y, z ) un sistema inercial y (u, v, w) un sistema m ovil que comparte origen con el primero. Los ejes z y w coinciden y el sistema (u, v, w) gira con velocidad angular constante alrededor del eje z . Entre ambos sistemas existe la relaci on geom etrica:

cos sen 0 x u cos 0 v y = sen 0 0 1 z w

(37)

Si consideramos el movimiento de un proyectil sometido s olo a la fuerza de la gravedad, en el sistema inercial sus ecuaciones del movimiento son: x = 0 y = 0 z = g (38)

19

Es una cuesti on simplemente algebraica escribir estas ecuaciones en funci on de las coordenadas del sistema m ovil: v 2 u = 0 u 2 u 2v = 0 v + 2 w = g (39)

Introduciendo variables de estado x1 = u, x2 = v , x3 = w, x4 = u , x5 = v y x6 = w : x 1 = x4 x 2 = x5 x 3 = x6 2 x1 + 2x 5 x 4 = 2 x2 2x 4 x 5 = x 6 = g (40)

Aparentemente, tenemos un sistema de seis ecuaciones, pero la tercera y la sexta no dependen de ninguna de las otras: se encuentran desacopladas y por tanto pueden resolverse de forma independiente la primera, segunda, cuarta y quinta por un lado, y la tercera y sexta por otro. Tercera y sexta tienen soluci on trivial: x6 (t) = z = x6 (0) gt = z (0) gt 1 1 2 (0)t gt2 x3 (t) = x3 (0) + x6 (0)t gt = z (0) + z 2 2 de donde obtenemos el tiempo de vuelo: t= 2z (0) g (42)

(41)

Integraremos durante este intervalo el sistema formado por las ecuaciones restantes. Como dijimos antes, si el sistema es lineal existe una soluci on anal tica y no es preciso hacer simulaci on alguna. Pero aqu no estamos en ning un momento intentando resolver los sistemas que vamos encontrando: nuestro inter es est a en el modelado y simulaci on.

20

Hemos simulado dos disparos, a lo largo del eje u del sistema de referencia m ovil (jo respecto a la plataforma, que gira), uno en sentido positivo y otro en sentido negativo. La Figura 9 representa la proyecci on de la trayectoria sobre el plano (u, v ). El resultado es que, con una velocidad inicial de 100 m/s y un angulo de lanzamiento de /4, para un alcance de unos 1000 m. se produce una desviaci on de aproximadamente 2.5 m. a la derecha, tomando como velocidad angular la mitad de la velocidad angular de rotaci on de la Tierra:
3

v(u)

3 1500

1000

500

500

1000

1500

Figura 9

9.

Conclusiones

Hemos usado la F sica elemental para plantear cinco versiones de un mismo sistema: un proyectil que se mueve bajo la acci on de la gravedad. No hemos resuelto anal ticamente ninguno de los problemas planteados, sino que hemos escrito algunos programas, implementando el algoritmo de Runge-Kutta de cuarto orden, para simularlos. Hemos prestado algo de atenci on al proceso de representaci on de los resultados, desde el archivo de datos en bruto generado por el programa de simulaci on hasta la inserci on de una gr aca en un documento .pdf (este documento). Se observar a que algunas veces la F sica es suciente para generar un modelo. Otras, ser a preciso hacer hip otesis razonables. En cualquier caso, hemos mostrado un m etodo susceptible de ser aplicado a

21

multitud de problemas din amicos, entendidos estos en el sentido amplio de sistemas descritos por variables que var an en el tiempo, no necesariamente sistemas mec anicos. Al mismo tiempo, esperamos haber sugerido otras preguntas, de las que dejamos constancia de unas cuantas: c omo modelar efectos aerodin amicos, como por ejemplo la sustentaci on del proyectil? qu e ocurre cuando se ha de tener en cuenta la curvatura de la Tierra? Qu e ocurre si un misil se aleja lo bastante de la supercie de la Tierra como para tener que considerar que la aceleraci on de la gravedad var a tanto en m odulo como en direcci on? Hay un ampl simo campo que ni siquiera hemos sugerido, aunque forma parte del curso: el control de sistemas. Relativas al control del sistema propuesto, algunas de las preguntas que podr amos hacer son: c omo conseguir que el proyectil impacte en una posici on predeterminada? c omo conseguir que siga una trayectoria prejada? c omo conseguir que el misil auto-corrija el efecto de perturbaciones atmosf ericas o errores en alguno de los par ametros de lanzamiento? En general c omo conseguir que la trayectoria del misil satisfaga alg un criterio prejado?

10.

Ap endice

En la p agina 9 se puede ver el c odigo del algoritmo de Runge-Kutta de cuarto orden aplicado al caso de un sistema de cuatro ecuaciones diferenciales ordinarias. A la vista de dicho c odigo, es f acil ver que la ampliaci on del algoritmo a un problema con un n umero arbitrariamente mayor de ecuaciones se reduce (desde el punto de vista de la codicaci on) a cortar y pegar. Esto motiva la pregunta de si ser a posible escribir una versi on de dicho algoritmo que tomase como par ametro el n umero de ecuaciones. Tambi en, si ser a posible aislar el algoritmo de las funciones concretas de los segundos miembros del sistema de ecuaciones diferenciales. La respuesta a ambas preguntas es s , y en este ap endice presentamos una implementaci on del algoritmo que toma como par ametros esencialmente el n umero de ecuaciones y un vector de punteros a funci on. Adem as, un puntero a un vector con las componentes del estado del sistema, los instantes inicial y nal, el intervalo de integraci on y un archivo donde se guardar an el tiempo y el estado del sistema en cada

22

instante. El algoritmo se presenta en la funci on rk4() y como puede verse en el programa de ejemplo que se adjunta el funcionamiento t pico consiste en: a) escribir las funciones de los segundos miembros (f0,1,2,3() en nuestro caso); b) Establecer condiciones iniciales; c) asignar los punteros a funci on y d) ejecutar rk4() /* Este programa ilustra el uso de la funcion rk4, que es un algoritmo generico de Runge-Kutta de cuarto orden que recibe como parametros el orden del sistema, un array de punteros a funcion, un vector de estado, instantes inicial y final y paso de integracion. Tambien un archivo, donde se escribiran en cada paso el valor del tiempo y las componentes del vector de estado. Los punteros a funcion apuntan a funciones que toman como parametros un puntero al vector de estado y quizas el tiempo. Seria mas modular escribir una funcion que calculase solo un paso de integracion, y usarla como base para una segunda que en un bucle efectuara la integracion completa. Pero como esta primera funcion tendria que reservar y liberar memoria en cada llamada, seria poco eficiente. Asi que escribimos en una unica funcion el bucle completo y luego, en cada caso, si es preciso anyadimos dentro de ese bucle lo que necesitemos. */

#include <stdio.h> #include <math.h> #include <malloc.h> #define pi 3.14159265 /* devuelve 1 si todo bien; 0 si no se pudo abrir el archivo donde se guardarn los datos */

23

int rk4(int n, double (*g[])(double*, double), double *x, double t0, double t1, double h, char *archivo) { int j; /* las ks: n filas por 4 columnas */ double *k=(double *)malloc(n*4*sizeof(double)); /* el estado siguiente */ double *y=(double *)malloc(n*sizeof(double)); /* la variable tiempo */ double t=t0; /* el archivo para los datos */ FILE *f; if ((f=fopen(archivo,"w"))==NULL) return(0); /* el bucle: aqui puede ponerse la condicion que se quiera */ while (t<t1-h){ /* paso cero: guardar tiempo y estado */ fprintf(f,"%10.4f ",t); for(j=0;j<n;++j) fprintf(f,"%10.4f ",*(x+j)); fprintf(f,"\n"); /* primer paso */ for(j=0;j<n;++j) *(k+4*j)=h*(*g[j])(x,t); /* segundo paso: preparar argumento y llamar */ for(j=0;j<n;++j) *(y+j)=(*(x+j)) + 0.5* (*(k+4*j)); for(j=0;j<n;++j) *(k+4*j+1)=h*(*g[j])(y,t+0.5*h); /* tercer paso: preparar argumento y llamar */ for(j=0;j<n;++j) *(y+j)=(*(x+j)) + 0.5* (*(k+4*j+1)); for(j=0;j<n;++j)

24

*(k+4*j+2)=h*(*g[j])(y,t+0.5*h); /* cuarto paso: preparar argumento y llamar */ for(j=0;j<n;++j) *(y+j)=(*(x+j)) + (*(k+4*j+2)); for(j=0;j<n;++j) *(k+4*j+3)=h*(*g[j])(y,t+h); /* quinto paso: calcular nuevo estado */ for(j=0;j<n;++j) *(x+j)=(*(x+j))+(1.0/6.0)*( (*(k+4*j)) + 2* (*(k+4*j+1)) + 2* (*(k+4*j+2)) + (*(k+4*j+3)) ); /* sexto paso: actualizar el tiempo */ t+=h; } /* cerrar el archivo */ fclose(f); /* liberar memoria */ free(k); free(y); return(1); } /* ejemplo de calculo de trayectorias de proyectiles pr.pdf, pag. 6 */ double f0(double *x, double t){ return(*(x+2)); } double f1(double *x, double t){ return(*(x+3)); } double f2(double *x, double t){ return(0); } double f3(double *x, double t){ return(-9.81);

25

} main() { int s; double (*g[4])(double *, double); /* valores iniciales */ double *x=(double *)malloc(4*sizeof(double)); *(x+0)=0; *(x+1)=0; *(x+2)=100*cos(pi/4); *(x+3)=100*sin(pi/4); g[0]=f0; g[1]=f1; g[2]=f2; g[3]=f3; s=rk4(4,g,x,0,10,0.1,"prueba.dat"); if (s) printf("\nOK\n"); else printf("\nNOK\n"); return; }

26