i

Introducci on a SCILAB
H ector Manuel Mora Escobar
Departamento de Matematicas Universidad Nacional de Colombia Bogota mayo del 2005

 INDICE GENERAL
1 INTRODUCCION 2 GENERALIDADES 2.1 Primeros pasos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Operaciones y funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Otros temas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Complejos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 3 3 6 9 9

2.5 Polinomios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.6 Lista de algunas herramientas . . . . . . . . . . . . . . . . . . . . . . . . 11 3 VECTORES Y MATRICES 13

3.1 Creacion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2 Notacion y operaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3 Funciones elementales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4 Solucion de sistemas de ecuaciones . . . . . . . . . . . . . . . . . . . . . 20 3.5 Otras funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4 PROGRAMAS 22

4.1 Guiones (scripts) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.2 Funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.3 Carpeta actual o por defecto . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.4 Comparaciones y operadores logicos . . . . . . . . . . . . . . . . . . . . . 27 4.5 Ordenes y control de ujo . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.5.1 if . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 ii

 INDICE GENERAL 4.5.2 4.5.3 4.5.4 4.5.5 4.6.1 4.6.2 4.6.3 4.6.4

iii

for . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 while . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 select . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Otras ordenes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Calculo numerico del gradiente . . . . . . . . . . . . . . . . . . . 33 Matriz escalonada reducida por las . . . . . . . . . . . . . . . . . 34 Aproximacion polinomial por m nimos cuadrados . . . . . . . . . 37 Factores primos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

4.6 Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.7 Miscelanea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5 GRAFICAS 43

5.1 Dos dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.2 Tres dimensiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.3 Otras funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.4 Creacion de un archivo Postscript . . . . . . . . . . . . . . . . . . . . . . 49 6 METODOS NUMERICOS 50

6.1 Sistemas de ecuaciones lineales . . . . . . . . . . . . . . . . . . . . . . . . 50 6.2 Solucion por m nimos cuadrados . . . . . . . . . . . . . . . . . . . . . . . 51 6.3 Una ecuacion no lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 6.4 Sistema de ecuaciones no lineales . . . . . . . . . . . . . . . . . . . . . . 53 6.5 Integracion numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.6 Derivacion numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.7 Interpolacion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 6.8 Aproximacion por m nimos cuadrados . . . . . . . . . . . . . . . . . . . . 56 6.9 Ecuaciones diferenciales ordinarias . . . . . . . . . . . . . . . . . . . . . . 58 6.10 Sistema de ecuaciones diferenciales . . . . . . . . . . . . . . . . . . . . . 58 7 OPTIMIZACION 7.1.1 60

7.1 Optimizacion lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Desigualdades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

iv 7.1.2 7.1.3

INDICE GENERAL Desigualdades y restricciones de caja . . . . . . . . . . . . . . . . 61 Igualdades, desigualdades y restricciones de caja . . . . . . . . . . 62

7.2 Optimizacion cuadratica . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

7.3 Optimizacion no lineal .18382 . . . 0 . Td . . . . . . .ados0 . . . Td0 . . . T271 . 63 9.23982d (y)Tj 1 0 j 0.23999 -3 T . 7.3 Optimizacion no

Cap tulo 1 INTRODUCCION

Este documento es una peque~ na introduccion a Scilab. Esta lejos de ser exhaustivo con respecto a los temas, es decir, muchos de los temas y posibilidades de Scilab no son tratados aqu . Ademas, tampoco es exhaustivo con respecto al contenido de cada tema, solo estan algunos aspectos de cada tema. El autor estara muy agradecido por los comentarios, sugerencias y correcciones enviados a: [email protected] o [email protected] Scilab fue desarrollado en el INRIA, Institut National de Recherche en Informatique et Automatique, un excelente instituto frances de investigacion. Posteriormente colaboro la escuela de ingenieros ENPC, Ecole Nationale de Ponts et Chaussees. Actualmente hay un gran consorcio con empresas como Renault, Peugeot-Citroen, CEA Commissariat a lEnergie Atomique, CNES Centre National dEtudes spatiales, Dassault Aviation, EDF Electricite de France, Thales... Sus principales caracter sticas son: software para calculo cient co interactivo programable de libre uso, con la condicion de siempre hacer referencia a sus autores disponible para diferentes plataformas: Windows, Linux, Sun, Alpha, ... El sitio o cial de Scilab es www.scilab.org 1

CAP ITULO 1. INTRODUCCION

All se encuentra informacion general, manuales, FAQs (frequent asked questions), referencias sobre reportes, diferencias con Matlab, lista de errores, ... All se puede \bajar" la version binaria o las fuentes para las diferentes plataformas. Un libro bastante completo sobre el tema es: Gomez C. ed., Engineering and Scientic Computing with Scilab, Birkhauser, Boston, 1999. Otros libros son: Allaire G y Kaber S.M., Introduction a ` Scilab, Exercises pratiques dalg` ebre lin eaire, Ellipses, Paris, 2002. Guerre-Delabriere S. y Postel M., M ethodes dapproximation : Equations di erentielles - Applications Scilab, niveau L3, Ellipses, Paris, 2004. Yger A., Th eorie et analyse du signal : Cours et initiation pratique via MATLAB et SCILAB (1 decembre 1999) Ellipses, Paris, 1999. Urroz G., Numerical and Statistical Methods With Scilab for Science and Engineering, vol. 1, 2, Booksurge, 2001. Gomez C., Bunks C. et al., Integrated Scientic Computing with SciLab, Birkhauser, Boston, 1998 Chancelier J.P. et al., Introduction a ` SCILAB (Collection IRIS) Springer, Paris, 2002 La utilizacion es igual para las diferentes plataformas, pero obviamente hay algunas diferencias intr nsecas al usar una u otra plataforma. Este peque~ no manual se re ere principalmente a la version 3.0 para Linux y para Windows (la ultima, a la fecha de hoy: 1 de abril del 2005). Deseo agradecer a todos los que tuvieron que ver con la creacion y actualizacion de Scilab. Realmente es una herramienta de uso libre muy poderosa. Me ha sido muy util para mi trabajo personal y para mis clases en la Universidad. Tambien deseo manifestar mi agradecimiento a todas las personas que de una u otra forma me han ayudado a corregir, modi car, actualizar y mejorar este documento, en particular, a Manuel Mej a.

Cap tulo 2 GENERALIDADES

Al activar Scilab aparece en la pantalla algo semejante a:

------------------------------------------Scilab-3.0 Copyright (c) 1989-2004 Consortium Scilab (INRIA, ENPC) -------------------------------------------

Startup execution: loading initial environment --> A partir de ese momento se puede escribir al frente de --> las ordenes de Scilab. Estas deben ser acabadas oprimiendo la tecla Enter , denominada tambien Intro o
| simplemente .

2.1

Primeros pasos

La orden
| --> t = 3.5

crea la variable t y le asigna el valor 3.5 . Ademas Scilab muestra el resultado de la orden desplegando en pantalla lo siguiente: 3

4 t = 3.5 -->

CAP ITULO 2. GENERALIDADES

De ahora en adelante, en este manual no se indicar a la tecla | ni tampoco el prompt --> . Por lo tanto deben sobreentenderse. Generalmente tampoco se mostrara el resultado desplegado en la pantalla por Scilab ante una orden recibida. Al dar la orden T = 4.5; se crea otra variable diferente con nombre T y se le asigna el valor 4.5. Scilab diferencia las letras minusculas de las mayusculas. La presencia de punto y coma al nal de la orden hace que Scilab no muestre el resultado en la pantalla. Sin embargo, la orden tuvo efecto. Scilab no hace diferencia entre numeros enteros y numeros reales. Los numeros se pueden escribir usando la notacion usual o la notacion cient ca. Por ejemplo, son validos los siguientes numeros. 3.5 -4.1234 3.14e-10 3.14E-10 0.0023e20 -12.345e+12 Al usar la orden who Scilab muestra las variables que esta usando en ese momento. Su respuesta es algo semejante a:
your variables are... T t startup %scicos_display_mode %scicos_short %helps percentlib soundlib s2flib roblib optlib autolib armalib alglib %F %T %z %t %f %eps ierr demolist scicos_pal MSDOS home xdesslib utillib metalib elemlib intlib mtlblib %s %nan %io %i

%scicos_menu PWD TMPDIR tdcslib siglib commlib polylib WSCI SCI %inf $ %e

2.1. PRIMEROS PASOS

using 5870 elements out of 1000000. and 48 variables out of

1791

Ah aparecen las variables t y T , de nidas anteriormente, y otras variables propias de Scilab. En las ordenes de Scilab los espacios en blanco antes y despues del signo igual no son indispensables. Se podr a simplemente escribir t=3.5 obteniendose el mismo resultado. Los espacios en blanco sirven simplemente para facilitar la lectura. Los nombres en Scilab constan hasta de 24 caracteres. El primero debe ser una letra o $ . Los otros pueden ser letras, numeros, #, _, $, ! . Por ejemplo: a12, pesoEsp, valor_ini Cuando en la orden no hay ninguna asignacion, sino simplemente una operacion valida, Scilab crea o actualiza una variable llamada ans . Por ejemplo, t+T produce el resultado ans 8. Las asignaciones pueden incluir operaciones con variables ya de nidas, por ejemplo x = t+T Si se da la orden t = 2*t se utiliza el antiguo valor de t para la operacion, pero, despues de la orden, t queda valiendo 7. Si se quiere conocer el valor de una variable ya de nida, basta con digitar el nombre de la variable y oprimir Enter. Es posible volver a repetir facilmente una orden dada anteriormente a Scilab. Utilizando las teclas correspondientes a las echas hacia arriba y hacia abajo, se obtiene una orden anterior y se activa oprimiendo | . Por ejemplo al repetir varias veces la orden x = (x+3/x)/2 se obtiene una aproximacion muy buena de 3. =

Tambien es posible, por medio de las echas (hacia arriba y hacia abajo), buscar una orden anterior para editarla y enseguida activarla.

CAP ITULO 2. GENERALIDADES

En una misma l nea de Scilab puede haber varias ordenes. Estas deben estar separadas por coma o por punto y coma. Por ejemplo, t1 = 2, t2 = 3; dt = t2-t1

2.2
+

Operaciones y funciones
* /

Los s mbolos sirven para las 4 operaciones aritmeticas. El signo - tambien sirve para indicar el inverso aditivo. Por ejemplo u = -t Para elevar a una potencia se utiliza el signo ^ u = 2^8, v = 2**8 ), como en la orden o tambien ** . Por ejemplo

Scilab utiliza para agrupar los parentesis redondos: ( x = (x+3/x)/2.

Puede haber parejas de parentesis metidas (anidadas) dentro de otras. En una expresion puede haber varios operadores. Las reglas de precedencia son semejantes a las de la escritura matematica usual. Los parentesis tienen prioridad sobre todos los operadores. Entre los operadores vistos hay tres grupos de prioridad. De mayor a menor, estos son los grupos de prioridad: ^ * + ** / -

Entre operadores de igual prioridad, se utiliza el orden de izquierda a derecha. Por ejemplo, 2*3+4^5-6/7 es equivalente a ((2*3)+(4^5))-(6/7). Scilab tiene prede nidas muchas funciones matematicas. A continuacion esta la lista de las funciones elementales mas comunes. abs : valor absoluto acos : arcocoseno acosh : arcocoseno hiperbolico asin : arcoseno asinh : arcoseno hiperbolico atan : arcotangente

2.2. OPERACIONES Y FUNCIONES atanh : arcotangente hiperbolica ceil : parte entera superior cos : coseno cosh : coseno hiperbolico cotg : cotangente coth : cotangente hiperbolica exp : funcion exponencial: ex fix : redondeo hacia cero (igual a int) floor : parte entera inferior int : redondeo hacia cero (igual a fix) log : logaritmo natural log10 : logaritmo decimal log2 : logaritmo en base dos max : maximo min : m nimo modulo : residuo entero rand : numero aleatorio round : redondeo sin : seno sinh : seno hiperbolico sqrt : raiz cuadrada tan : tangente tanh : tangente hiperbolica

El signi cado de la mayor a de estas funciones es absolutamente claro. La siguiente tabla muestra varios ejemplos utilizando las funciones de parte entera y redondeo. x ceil(x) floor(x) int(x) round(x) 2.0 2. 2. 2. 2. 1.8 2. 1. 1. 2. 1.5 2. 1. 1. 2. 1.2 2. 1. 1. 1. - 3.1 - 3. - 4. - 3. - 3. - 3.5 - 3. - 4. - 3. - 4. - 3.8 - 3. - 4. - 3. - 4.

CAP ITULO 2. GENERALIDADES

Otra funcion matematica, esta ya con dos parametros de entrada, es modulo. Sus dos parametros deben ser enteros. El resultado es el residuo de la division entera. modulo(17,5) da como resultado 2. Para tener informacion mas detallada sobre alguna funcion basta con digitar help y a continuacion el nombre de la funcion o de la orden. Por ejemplo help floor Obviamente se requiere que la funcion floor exista. Si no se conoce el nombre de la funcion, pero se desea buscar sobre un tema, se debe utilizar apropos. Por ejemplo: apropos polynomial da informacion sobre las funciones que tienen que ver con polinomios. En cambio, help polynomial informa que no hay manual para polynomial. Ademas de estas funciones elementales, Scilab tiene muchas mas funciones como las funciones de Bessel, la funcion gama, ... Mediante la barra de menu, con la opcion Help seguida de Help Dialog se obtiene un catalogo resumido de las herramientas de Scilab. La lista de funciones elementales de Scilab es la siguiente:

abs, acos, acosh, acoshm, acosm, addf, adj2sp, amell, and, asinh, asinhm, asinm, atan, atanh, atanhm, atanm, besseli, besselj, besselk, bessely, binomial, bloc2exp, bloc2ss, calerf, ceil, cmb_lin, conj, cos, cosh, coshm, cosm, cotg, coth, cothm, cumprod, cumsum, delip, diag, dlgamma, double, erf, erfc, erfcx, eval, eye, fix, floor, frexp, full, gamma, gammaln, gsort, imag, int, int16, int32, int8, integrate, interp, interpln, intersect, intsplin, inttrap, isdef, isinf, isnan, isreal, kron, ldivf, lex_sort, linspace, log, log10, log2, logm, logspace, max, maxi, mean, median, min, mini, minus, modulo, mps2linpro, mtlb_sparse, mulf, nnz, norm, not, ones, or, pen2ea, pertrans, pmodulo, prod, rand, rat, rdivf, real, round, sign, signm, sin, sinh, sinhm, sinm, size, smooth, solve, sort, sp2adj, sparse, spcompack, speye, spget, splin, spones, sprand, spzeros, sqrt, sqrtm, squarewave, ssprint, ssrand, st_deviation, subf, sum, sysconv, sysdiag, syslin, tan, tanh, tanhm, tanm, toeplitz, trfmod, trianfml, tril, trisolve, triu, typeof, uint16, uint32, uint8, union, unique, zeros.

2.3. OTROS TEMAS

2.3

Otros temas

Se puede modi car el formato utilizado por Scilab para mostrar los resultados, mediante format. Si se da la orden format(16) a partir de ese momento, Scilab utilizara 16 \columnas" (16 posiciones) para mostrar cada numero. Estas 16 columnas incluyen el espacio para el signo la parte entera y el punto. Por defecto, Scilab usa 10 posiciones. format(e,14) La orden anterior sirve para utilizar notacion cient ca con 14 posiciones. Tambien se puede utilizar simplemente format(e) Para regresar al formato inicial, el formato \variable" (el prede nido por Scilab) se usa format(v) o, por ejemplo, format(v, 10) Scilab tiene prede nidas algunas constantes especiales cuyos nombres estan precedidos del signo % . Para los valores e, , 1, sus nombres son: %e , %pi , %i . Observe que dos de estas variables aparecieron al utilizar who. Despues de la siguiente asignacion, la variable r tendra el valor -1. r = log(1/%e)

2.4

Complejos

Scilab maneja de manera sencilla los numeros complejos. Estos pueden ser de nidos de varias maneras. Por ejemplo, suponiendo que r vale 5, las dos ordenes siguientes a = 3 + 4*r*%i b = sqrt(-4) Para las operaciones con numeros complejos (suma, resta, multiplicacion, ...) utilizan exactamente los mismos s mbolos + - * / ** ^. de nen dos variables, a = 3 20i, b = 2i. se

Las funciones real, imag y conj permiten obtener la parte real, la parte imaginaria y el conjugado de un complejo. Si se utiliza la funcion abs con un complejo, se obtiene la magnitud o modulo de el. Las funciones de Scilab usadas para funciones reales elementales que tienen generalizaciones en complejos, se pueden usar tambien para los complejos, por ejemplo, sin,

10 cos, log, ... As , es completamente l cito z = 3 + 4*%i; r = sin(z) El resultado mostrado por Scilab en pantalla es r = 3.853738 - 27.016813i

CAP ITULO 2. GENERALIDADES

2.5

Polinomios

Un polinomio se puede de nir de dos maneras: por sus coe cientes o por sus ra ces. Es necesario ademas indicar la variable simbolica para el polinomio. La orden p = poly([2 3 5 7], "x", "coeff") de ne en la variable p el polinomio 2 + 3x + 5x2 + 7x3 . La orden q = poly([2 3 5], "x", "roots") de ne en la variable q el polinomio 30+31x 10x2 + x3 cuyas ra ces son exactamente 2, 3 y 5. Escribir q = poly([2 3 5], "x") produce exactamente el mismo resultado, o sea, "roots" es el tipo de de nicion por defecto. La doble comilla " puede ser reemplazada por la comilla sencilla . Mas aun, se puede reemplazar coeff por c y roots por r . Es l cito escribir r = poly([6 7 8], y, c) . La funcion roots calcula las ra ces de un polinomio, sean estas reales o complejas. Por ejemplo roots(p) Con polinomios se pueden hacer sumas, multiplicaciones, restas, multiplicacion por un numero. Deben ser polinomios en la misma variable. Por ejemplo: v = p + q + p*q - 3.1*q Tambien se puede elevar un polinomio a una potencia, por ejemplo, r = p^3 La funcion coeff tiene dos parametros, el primero es el polinomio y el segundo la potencia. La siguiente orden asigna a la variable k el valor 10, el coe ciente de x 2 en el polinomio q. k = coeff(q, 2) Si se utiliza simplemente c = coeff(q)

2.6. LISTA DE ALGUNAS HERRAMIENTAS

11

se obtendran todos los coe cientes. La variable c sera un vector (ver cap tulo siguiente sobre matrices y vectores). Si se utiliza p = poly(a,x), donde a es una matriz cuadrada (ver cap tulo siguiente sobre matrices y vectores), se obtiene el polinomio caracter stico de de la matriz a. Para evaluar un polinomio p en un valor t se usa horner(p, t) Por ejemplo horner(q, 1) dara como resultado 8. Si q es un polinomio, es l cito utilizar la orden r = horner(p, q) para obtener p(q(x)).

2.6

Lista de algunas herramientas

Scilab tiene numerosas herramientas para diferentes temas. A continuacion hay una lista de ellas. Algebra lineal (cap tulo 3) Gra cas (cap tulo 5) Optimizacion METANET: grafos y redes Analisis y control de sistemas lineales Procesamiento de se~ nales Modelacion y simulacion ARMA (autoregressive moving average) Simulacion SCICOS: modelamiento y simulacion de sistemas dinamicos h bridos Funciones de entrada y salida Manejo de funciones y librer as Manipulacion de cadenas de caracteres Dialogos: ventanas y botones Calculos con polinomios

12 Funciones de distribucion Control robusto PVM: parallel virtual machine Traducciones de lenguajes y datos GeCI: comunicacion con otras aplicaciones Interfaz con Maple

CAP ITULO 2. GENERALIDADES

Cap tulo 3 VECTORES Y MATRICES

En Scilab no hay vectores como tales, los vectores se deben asimilar a matrices de una sola la o vectores la (tama~ no 1 n) o a matrices de una sola columna o vectores columna (tama~ no n 1).

3.1

Creaci on
11 12 13 14 15 21 22 23 24 25 31 32 33 34 35

La matriz

se puede de nir por medio de o tambien por medio de

a = [ 11 12 13 14 15; 21 22 23 24 25; 31 32 33 34 35] a = [ 11,12,13,14,15; 21,22,23,24,25; 31,32,33,34,35] o por a = [ 11 12 13 14 15 21 22 23 24 25 31 32 33 34 35] Scilab mostrara el siguiente resultado en la pantalla: a ! ! ! = 11. 21. 31. 12. 22. 32. 13. 23. 33. 14 24 34 15. ! 25. ! 35. ! 13

14

CAP ITULO 3. VECTORES Y MATRICES

La de nicion de la matriz se hace por las. Los elementos de una misma la se separan por medio de espacios en blanco o por medio de comas. Una la se separa de la siguiente por medio de punto y coma o por medio de cambio de l nea. Scilab permite crear rapidamente algunos tipos especiales de matrices: ones(4,5) es una matriz de unos de tama~ no 4 5 zeros(4,5) es una matriz de ceros de tama~ no 4 5 eye(4,4) es la matriz identidad de orden 4

rand(20,30) es una matriz aleatoria de tama~ no 20 30

Algunas veces es util, especialmente para gra cas de funciones (tema que se vera en un cap tulo posterior), crear vectores con elementos igualmente espaciados, por ejemplo x = [1 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6] Esto tambien se puede hacer mas facilmente as : x = 1:0.2:2.6 Los valores 1.0, 0.2 y 2.6 corresponden respectivamente al l mite inferior, al incremento y al l mite superior. Si no se especi ca incremento se supone que es uno. Escribir y = 2:9 es equivalente a y = 2:1:9. Cuando se desea una columna basta con utilizar el operador de trasposicion, por ejemplo, z = (1:0.2:2.6)

3.2

Notaci on y operaciones

a(2,5) denota el elemento o entrada de a en la la 2 y en la columna 5. a(3,:) denota la tercera la de la matriz a. a(:,4) denota la cuarta columna de la matriz a. a(1:2,2:5) denota la submatriz de tama~ no 2 4, formado por los elementos que estan en las las 1, 2 y en las columnas 2, 3, 4, 5, o sea, la matriz 12 13 14 15 22 23 24 25 La anterior de nicion de submatrices se puede generalizar utilizando vectores( la o columna). Por medio de las siguientes ordenes u = [ 1 2 1], v =[2 4], aa = a(u,v) 12 14 22 24 12 14 se asigna a la variable aa la matriz

 Y OPERACIONES 3.2. NOTACION

15

Si a es una matriz m n, entonces a(:) es un vector columna de tama~ no mn 1, obtenido colocando primero la primera columna de a, enseguida la segunda columna, luego la tercera y as sucesivamente. Si x es un vector columna, entonces x(:) es el mismo vector columna. Si x es un vector la, entonces x(:) es lo mismo que x. Dicho de otra forma, si x es un vector ( la o columna), entonces x(:) es un vector columna con las mismas entradas de x. As , x = x(:) convierte el vector x en un vector columna. Si x es un vector la, se puede escribir x(3) en lugar de x(1,3). Analogamente, si y es un vector columna, se puede escribir y(2) en lugar de y(2,1). Por medio de * se puede hacer el producto entre un numero y una matriz. Los signos + y - permiten hacer sumas y restas entre matrices del mismo tama~ no. Cuando el producto de matrices es posible (numero de columnas de la primera matriz igual al numero de las de la segunda), este se indica por medio de * . La transposicion de una matriz se indica por medio de comilla, por ejemplo, c = a*a crea la matriz c , producto de la traspuesta de a y de a . Por construccion c es una matriz cuadrada. La mayor a de las funciones de Scilab utilizadas para numeros reales se pueden aplicar a matrices. En este caso se obtiene una matriz del mismo tama~ no, en la que se ha aplicado la funcion a cada uno de los elementos de la matriz. Por ejemplo, las dos ordenes siguientes son equivalentes. D = sin([1 2; 3 4]) , D = [sin(1) sin(2); sin(3) sin(4)] Para matrices (y vectores) del mismo tama~ no se puede hacer la multiplicacion elemento a elemento utilizando .* . De manera analoga se puede hacer la division elemento a elemento. Por ejemplo: P = rand(3,5), H = a.^3 Si a es una matriz rectangular a.^(1/2) es lo mismo que sqrt(a) es decir, una matriz del mismo tama~ no, cuyas entradas son las raices cuadradas de las entradas de a. En cambio, si a es una matriz cuadrada, a1 = a^(1/2) Q = rand(3,5), r = P.*Q Tambien es posible elevar a una potencia los elementos de una matriz, por ejemplo,

16

CAP ITULO 3. VECTORES Y MATRICES

es una matriz raiz cuadrada de a, o sea, a1*a1 es lo mismo que a. Como se vera mas adelante, para gra car se requieren dos vectores, uno con los valores de la variable x y otro (del mismo tama~ no) con los valores de la variable y . Las dos ordenes siguientes permiten crear los dos vectores para un polinomio: x = (-2:0.01:3); y = 3*x.^4 + x.^3 - 5*x.*x + 3.14; Para matrices cuadradas, es posible calcular directamente una potencia. Por ejemplo, G = rand(6,6)^3 Si los tama~ nos son compatibles, dos o mas matrices se pueden \pegar" para obtener una matriz de mayor tama~ no, por ejemplo, AA = [a rand[3,2)] B = [rand(2,5); a; eye(5,5)] Como a fue de nida de tama~ no 3 5, entonces AA es de tama~ no 3 7 y B es de tama~ no 10 5. Si x es un vector, la orden x(2) = [] suprime la segunda componente y desplaza el resto. Por ejemplo, x = [2 3 5 7 11]; x(2) = [] produce el resultado x ! = 2. 5. 7. 11. ! La orden y = [] permite crear la matriz y de tama~ no 0 0.

De manera analoga, si a es un matriz, la orden a(:,2) = [] suprime la segunda columna.

3.3

Funciones elementales

Hay numerosas funciones de Scilab para el manejo de matrices. Algunas de las mas usadas son: rank(a) calcula el rango de a. det(c) determinante de una matriz cuadrada c. inv(c) inversa de una matriz cuadrada e invertible c.

3.3. FUNCIONES ELEMENTALES rref(a) matriz escalonada reducida por las equivalente a a. expm(C) produce la matriz eC , donde C es una matriz cuadrada. 1 1 1 eC = I + C + C 2 + C 3 + C 4 + 2 6 24

17

diag(c) produce un vector columna con los elementos diagonales de la matriz cuadrada c. diag(x) produce una matriz diagonal con los elementos del vector ( la o columna) x. y = sort(x) ordena el vector x de manera decreciente. [y, k] = sort(x): y es el vector ordenado de manera decreciente, k es un vector que contiene los ndices del ordenamiento, o sea, y = x(k). b = sort(a) ordena la matriz a de manera decreciente, considerando cada matriz como un vector formado por la primera columna, la segunda columna, ..., la ultima columna. a= 3.2 3.1 , 3.4 3.8 b= 3.8 3.2 3.4 3.1

b = sort(a, r) ordena la matriz a de manera decreciente por columnas. Atenci on, aunque r tiene un signi cado interno de las, el resultado externo es un ordenamiento de las columnas. b= 3.4 3.8 3.2 3.1

b = sort(a, c) ordena la matriz a de manera decreciente por las. Atenci on, aunque c tiene un signi cado interno de columnas, el resultado externo es un ordenamiento de las las. b= 3.2 3.1 3.8 3.4

m = max(x) calcula el maximo del vector ( la o columna) x. [m, k] = max(x): esta el maximo. m es el maximo del vector x, k indica la posicion donde

m = max(a) calcula el maximo de la matriz a.

18

CAP ITULO 3. VECTORES Y MATRICES [m, k] = max(a): m es el maximo de la matriz a, k es un vector 1 2 e indica la posicion donde esta el maximo. m = max(a, r) : columnas de a. m es un vector la (row) que contiene los maximos de las

[m, k] = max(a, r) : m es un vector la que contiene los maximos de las columnas de a, k es un vector la que contiene las posiciones de los maximos. min semejante a max pero para el m nimo. m = mean(x) calcula el promedio del vector ( la o columna) x. m = mean(a) calcula el promedio de la matriz a. m = mean(a, r) : columnas de a. m es un vector la (row) que contiene los promedios las

m = mean(a, c) : m es un vector columna que contiene los promedios las las de a. median semejante a mean pero para la mediana. st-deviation semejante a mean pero para la desviacion estandar. sum semejante a mean pero para la suma. prod semejante a mean pero para el producto. norm(x) calcula la norma euclidiana del vector x ( la o columna). calcula la norma lp del vector x:
n 1/p

norm(x, p)

i=1

|xi |

norm(x, inf)

calcula la norma \del maximo" del vector x:

1in

max |xi |.

norm(a) = norm(a, 2) calcula, para la matriz a, la norma matricial generada por la norma euclidiana, o sea, el mayor valor singular de a. norm(a, 1) calcula, para la matriz a, la norma matricial generada por la norma l1 , o sea,
m 1j n

max

i=1

|aij |.

3.3. FUNCIONES ELEMENTALES

19

norm(a, inf) calcula, para la matriz a, la norma matricial generada por la norma l , o sea,
n 1im

max

j =1

|aij |.

norm(a, fro) calcula, para la matriz a, la norma de Frobenius, o sea,

m n 1/2

a2 ij
i=1 j =1

La orden fc = size(a) proporciona un vector la 1 2 con el numero de las y el numero de columnas de a. La orden size(a,1) proporciona el numero de las de a. Analogamente, size(a,2) proporciona el numero de columnas de a. Por medio de tril se puede obtener la parte triangular inferior (low) de una matriz cuadrada. Por ejemplo, a = [1 2 3; 4 5 6; 7 8 9]; L = tril(a) produce el resultado L ! ! ! = 1. 4. 7. 0. 5. 8. 0. ! 0. ! 9. !

La funcion triu produce la parte triangular superior (upper). Si a es una matriz cuadrada, por medio de v = spec(a) se obtiene un vector columna con los valores propios (reales o complejos) de a. Con la orden [v, s] = spec(a) se obtiene la matriz v cuyas columnas son vectores propios y la matriz diagonal s cuyos elementos diagonales son los valores propios. Si A es una matriz simetrica y de nida positiva, entonces se puede factorizar en la forma A = U T U , donde U es una matriz triangular superior de diagonal positiva. Esta es la llamada factorizacion de Cholesky. Esta matriz se puede obtener por U = chol(A)

20

CAP ITULO 3. VECTORES Y MATRICES La factorizacion QR de una matriz A es la expresion A = QR donde Q es una matriz ortogonal y R es una matriz del tama~ no de A, \triangular " (no necesariamente cuadrada). Por ejemplo, [q, r] = qr(a) La factorizacion LU de una matriz cuadrada invertible A es la expresion P A = LU donde P es una matriz de permutacion, L es una matriz triangular inferior con unos en la diagonal y U es una matriz triangular superior. Por ejemplo, [L, U, P] = lu(A) La descomposicion SVD, descomposicion en valores singulares, de una matriz A es la expresion A = U DV T donde U y V son matrices ortogonales y D es una matriz del

ma2 0 Td (or338 Aonej cua Sol 0

3.5. OTRAS FUNCIONES

21

Teoricamente, el resultado debe ser el mismo. Desde el punto de vista de precision numerica, los resultados son semejantes. Para matrices medianamente grandes hay diferencias en tiempo. La primera forma gasta, aproximadamente, tres veces mas tiempo que la segunda. Fuera del caso de solucion unica hay otros dos casos: el caso de sistema inconsistente (sin solucion) y el caso de muchas soluciones (numero in nito de soluciones). Si el sistema es inconsistente, entonces se desea encontrar una seudosolucion, la solucion por m nimos cuadrados, o sea, se busca x que minimice ||Ax b||2 2 . En este caso la orden x = A\b encuentra una de estas\soluciones". Si las columnas de A son linealmente independientes, esta seudosolucion es unica. Para mayor informacion, use la ayuda: help backslash

3.5

Otras funciones

La lista de funciones de Scilab para Algebra Lineal es la siguiente: aff2ab, balanc, bdiag, chfact, chol, chsolve, classmarkov, coff, colcomp, companion, cond, det, eigenmarkov, ereduc, exp, expm, fstair, fullrf, fullrfk, genmarkov, givens, glever, gschur, gspec, hess, householder, im_inv, inv, kernel, kroneck, linsolve, lu, ludel, lufact, luget, lusolve, lyap, nlev, orth, pbig, pencan, penlaur, pinv, polar, proj, projspec, psmall, qr, quaskro, randpencil, range, rank, rcond, rowcomp, rowshuff, rref, schur, spaninter, spanplus, spantwo, spchol, spec, sqroot, sva, svd, sylv, trace

Cap tulo 4 PROGRAMAS

En Scilab hay dos tipos de programas: los guiones o libretos (scripts) y las funciones. Un guion es simplemente una secuencia de ordenes de Scilab. No tiene parametros (\argumentos") de entrada ni de salida. En cambio una funcion s los tiene. Por otro lado, las variables de nidas en un guion son globales, es decir, despues del llamado del guion estas variables siguen existiendo. En cambio en una funcion, las variables de nidas dentro de la funcion dejan de existir una vez nalizada la ejecucion de la funcion, son variables locales.

4.1

Guiones (scripts)

Un guion es simplemente un archivo ASCII en el que hay una sucesion de ordenes de Scilab. Generalmente tiene la extension .sce pero eso no es obligatorio. Puede estar colocado en cualquier carpeta. En el ejemplo que sigue se va a hacer lo siguiente: crear aleatoriamente una matriz, crear aleatoriamente la solucion x0 , crear los terminos independientes correspondientes a x0 , hallar la solucion, calcular la norma del error cometido. Scilab viene con un editor llamado Scipad. Se activa mediante la opcion Editor de la barra de menu de Scilab. Ya sabiendo lo que se va a hacer, con Scipad o con un editor cualquiera, creamos el archivo c:\misci\ensayo01.sce 22

4.1. GUIONES (SCRIPTS) cuyo contenido sea el siguiente: n = 100; A = rand(n,n); x0 = rand(n,1); b = A*x0; x = A\b; residuo = norm(x-x0)

23

Obviamente el nombre de la carpeta misci es simplemente un ejemplo y hubiera podido ser cualquier otra carpeta. Una vez guardado el contenido, desde el ambiente Scilab se da la orden exec c:\misci\ensayo01.sce Esto hace que se ejecuten todas las ordenes contenidas en el archivo. Mediante who, o de cualquier otra forma, se puede veri car que las variables n, A, x0, ... fueron creadas y todav a existen. Dar la orden exec c:\misci\ensayo01.sce tambien se puede hacer en Windows por medio de la barra de menu con las opciones File y despues Exec. Subiendo y bajando de nivel se busca la carpeta adecuada y se hace doble clic con el boton derecho del raton en el archivo ensayo01.sce. En Linux se hace de manera muy parecida: en la barra de menu File, enseguida File Operations y nalmente Exec despues de haber escogido el archivo deseado. Si se usa SciPad, se obtiene el mismo resultado con la barra de menu (de SciPad) mediante la opcion Execute y despues Load into Scilab. Si se desea, se puede editar el archivo, hacer algunos cambios, guardarlos y de nuevo activar el guion mediante exec ... En el siguiente ejemplo suponemos que se tiene una matriz A de tama~ no m n, m n, r (A) = n, es decir, tiene mas las que columnas y sus columnas son linealmente independientes. Consideremos el siguiente sistema de ecuaciones A x = b. Probablemente este sistema de ecuaciones no tiene solucion en el sentido estricto. Su seudosolucion o solucion por m nimos cuadrados esta dada por la solucion de la ecuacion normal AT A x = AT b. proveniente de minimizar el cuadrado de la norma del error: = ||A x b||2 2. Sea c:\estad\ensayo02.sce el archivo que de ne los datos y lleva a cabo este proceso. Su contenido puede ser:

24 // solucion por minimos cuadrados // a = [ 1 2; 3 4; 5 6]; b = [ 3 7 10]; // x = (a*a)\(a*b) e = norm(a*x-b)^2

CAP ITULO 4. PROGRAMAS

Las l neas que empiezan por // son l neas de comentarios (como en C++). Obviamente para activar este otro archivo se necesita dar la orden exec c:\estad\ensayo02.sce De manera natural aparece la pregunta: C omo hacer el mismo proceso con otro sistema de ecuaciones sin tener que cambiar el archivo? Las funciones son la respuesta.

4.2

Funciones

En un archivo ASCII puede haber varias funciones. Generalmente el nombre de los archivos de funciones tienen la extension .sci. El esquema general de una funcion es el siguiente: function [res1, res2, ...] = nombrefuncion(par1, par2, ...) ... ... endfunction El signi cado de res1 es resultado 1 o tambien parametro de salida 1. El signi cado de par2 es parametro de entrada 2 o simplemente parametro 2. Cuando la funcion tiene un unico resultado, el esquema puede ser simplemente: function result = nombrefuncion(par1, par2, ...) ... ... endfunction El siguiente archivo llamado c:\misci\misfunc.sci tiene varias funciones function [x, error2] = seudoSol(A, b) // solucion por minimos cuadrados x = (A*A)\(A*b) error2 = norm(A*x-b)^2

4.2. FUNCIONES endfunction //-------------------------------------------------------function [x, y] = polarCart(r, t) // Conversion de coordenadas polares a cartesianas. x = r*cos(t) y = r*sin(t) endfunction //-------------------------------------------------------function [x, y] = polarCartGr(r, t) // Conversion de coordenadas polares a cartesianas, // el angulo esta dado en grados. [x, y] = polarCart(r, t*%pi/180) endfunction //-------------------------------------------------------function fx = f(x) fx = x(1)^2 + x(2)^2 endfunction Este archivo de funciones se debe cargar mediante la orden getf c:\misci\misfunc.sci

25

o tambien, en Windows, mediante las opciones de la barra de menu File y despues Exec. En Linux, en la barra de menu File , File Operations y nalmente Getf despues de haber escogido el archivo deseado. Tambien con la barra de menu de SciPad: Execute y Load into Scilab . Una vez cargado el archivo, las funciones se pueden utilizar como las otras funciones de Scilab. Por ejemplo, son validas las siguientes ordenes: [x1, y1] = polarCart(2, 0.7854) [u, v] = polarCartGr(3, 30) a = [1 2; 3 4; 5 6]; b = [3 7 10]; [x,err] = seudoSol(a,b) [x, dif] = seudoSol([1 2; 3 4; 5 6], [3 7 10]) valor = f([3; 4]) x = [5; 6], res = f(x) Cuando una funcion produce mas de un resultado, tambien se puede utilizar asignando menos de los resultados previstos. Por ejemplo, la utilizacion completa de polarCartGr puede ser [a, b] = polarCartGr(3,30) lo cual produce el resultado b =

26

CAP ITULO 4. PROGRAMAS

1.5 = 2.5980762

En cambio, la orden c = polarCartGr(3,30) produce el resultado c = 2.5980762 Otra caracter stica de las funciones es que una funcion puede llamar una o varias funciones. Observese que la funcion polarCartGr utiliza la funcion polarCart. En los ejemplos anteriores de funciones, hay simplemente una secuencia de ordenes que siempre se repite de la misma manera. Con mucha frecuencia esto no es as . Normalmente dentro de una funcion, dependiendo de los datos o de resultados intermedios, hay que tomar algunas decisiones, repetir un proceso, abortar un proceso, ... Esto se logra mediante los operadores relacionales, los operadores logicos y las estructuras de control. Esto se vera en secciones posteriores de este mismo cap tulo.

4.3

Carpeta actual o por defecto

Al arrancar, Scilab tiene de nido un subdirectorio (o carpeta) preferencial, actual o por defecto. Para saber el nombre de esta carpeta, se da la orden pwd La respuesta de Scilab puede ser: ans =

C:\WINDOWS\Escritorio En Linux la respuesta podr a ser ans =

/home/h

 4.4. COMPARACIONES Y OPERADORES LOGICOS

27

Esto quiere decir que si un archivo de funciones o un archivo tipo guion estan ubicados all , no es necesario, al utilizar getf o exec, escribir la ubicacion completa del archivo, basta dar el nombre del archivo. Por ejemplo, si el archivo ejemplo4.sce esta en la carpeta C:\WINDOWS\Escritorio, no es necesario, dar la orden exec \WINDOWS\Escritorio\ejemplo4.sce Basta digitar exec ejemplo4.sce Para decir a Scilab que cambie la carpeta por defecto, se usa la orden chdir, por ejemplo, chdir c:\algebra\matrices\

4.4

Comparaciones y operadores l ogicos

< <= > >= == ~= <> & | ~ menor menor o igual mayor mayor o igual igual diferente diferente y o no

4.5

Ordenes y control de ujo

La principales estructuras de control de Scilab son: if select y case for while Otras ordenes relacionadas con el control de ujo son: break

28 return equivalente a resume

CAP ITULO 4. PROGRAMAS

En este manual introductorio se supone que el lector tiene alg un conocimiento de un lenguaje de programaci on. Aqu no se explica como es el funcionamiento de la estructura while. Simplemente hay algunos ejemplos sobre su escritura en Scilab.

abort

4.5.1

if

Una forma sencilla de la escritura de la escritura if es la siguiente: if condici on then ... ... end La palabra then puede ser reemplazada por una coma o por un cambio de l nea. Entonces se puede escribir if condici on ... ... end Obviamente tambien existe la posibilidad else: if condici on ... else ... end Estas estructuras se pueden utilizar directamente dentro del ambiente interactivo de Scilab. Por ejemplo, se puede escribir directamente, sin necesidad de un guion if x < 0, fx = x*x, else fx = x*x*x, end Resulta ine ciente hacer un guion unicamente para las ordenes de la l nea anterior. Pero, cuando hay muchas ordenes y muchos controles y es necesario depurar el proceso, es casi indispensable hacer un guion o una funcion.

 4.5. ORDENES Y CONTROL DE FLUJO

29

4.5.2

for

La estructura for tiene la siguiente forma: for var = lim1:incr:lim2 ... ... end Esto quiere decir que la variable var empieza en el l mite inferior lim1, despues va incrementar el valor en el incremento incr (puede ser negativo). Si el incremento es 1, este se puede suprimir. A continuacion algunos ejemplos. for i = 2:3:14 ... ... end for j = 2:-3:-10 ... ... end for k = 2:8 ... ... end Una estructura for puede estar anidada dentro de otro for o dentro de un if.

4.5.3

while

La forma general es: while condici on ... ... end Por ejemplo e = 1; while e+1 > 1 e = e/2; end

30

CAP ITULO 4. PROGRAMAS

4.5.4

select

La forma general es: select variable case valor1 then ... ... case valor2 then ... ... case valor3 then ... ... case valor4 then ... ... ... else ... ... end La palabra then puede ser reemplazada por una coma o por un cambio de l nea. La parte else es opcional. Entonces se puede escribir select variable case valor1 ... ... case valor2, ... case valor3 ... ... case valor4

 4.5. ORDENES Y CONTROL DE FLUJO ... ... end Por ejemplo select indic case 0, a = b; case 1 a = b*b; b = 1; case 3 a = max(a,b); b = a*a; else a = 0; b = 0; end

31

4.5.5

Otras o rdenes

La orden break permite la salida forzada (en cualquier parte interna) de un bucle for o de un bucle while. La orden return permite salir de una funcion antes de llegar a la ultima orden. Todos los par ametros de salida o resultados deben estar denidos con anterioridad. Otra orden que sirve para interrumpir una funcion, en este caso interrumpiendo la evaluacion, es abort. En la siguiente funcion se calcula el maximo comun divisor por el algoritmo de Euclides. function [maxcd, indic] = mcd(a, b) // Maximo comun divisor de a, b enteros positivos. // indic valdra 1 si se calculo el m.c.d // 0 si los datos son inadecuados. indic = 0 maxcd = 0 if round(a) ~= a | round(b) ~= b return end if a < 1 | b < 1 return

32 end if a < b t = a a = b b = t end indic = 1 while 1 == 1 r = modulo(a, b) if r == 0 maxcd = b return end a = b b = r end endfunction

CAP ITULO 4. PROGRAMAS

En el ejemplo anterior, el ultimo return esta anidado en un solo bucle y despues acaba la funcion, luego se podr a cambiar por un break. La condicion del while siempre es cierta, luego la salida del bucle se obtiene siempre en la mitad del cuerpo del while. En una funcion no es necesario el punto y coma despues de una orden, pues de todas maneras Scilab no muestra el resultado de la asignacion (en Matlab si es necesario el punto y coma). Si de nitivamente, en una funcion, se desea mostrar en pantalla algunos valores intermedios se debe hacer por medio de disp o por medio de printf. Por ejemplo, despues del calculo de r se puede utilizar la orden disp(a, b, r) Observe que primero escribe el valor de r, despues el de b y nalmente el de a. La orden printf es una emulacion de la orden del mismo nombre del lenguaje C. Por ejemplo se podr a utilizar la orden printf( a = %3d b = %3d r = %3d\n, a, b, r) Esta emulacion utiliza una comilla en lugar de la comilla doble. El s mbolo \n se utiliza para cambio de l nea. Otras ordenes que se pueden utilizar para escritura son: fprintf print write

4.6. EJEMPLOS

33

4.6

Ejemplos

Los siguientes ejemplos, son simplemente casos que pretenden ilustrar la manera de programar en Scilab. En ellos se busca mas la claridad que la e ciencia numerica.

4.6.1

C alculo num erico del gradiente

En este ejemplo se desea aproximar numericamente el gradiente de una funcion, que va de Rn en R, cuyo nombre es exactamente f. Cuando se desea cambiar de funcion, es necesario editar la funcion f e introducir los cambios necesarios. La aproximacion utilizada es f f (x + hei ) f (x hei ) , (x) xi 2h donde e1 , e2 , ..., en forman la base canonica. A continuacion aparece la funcion f y la funcion que aproxima el gradiente. //-------------------------------------------------------function fx = f(x) fx = 3*x(1)^2 + 5*x(2) endfunction //-------------------------------------------------------function g = gradf(x) // gradiente de la funcion f en el punto x //******** h = 0.001 //******** h2 = 2*h n = max(size(x)) g = zeros(n,1) for i =1:n xi = x(i) x(i) = xi+h f1 = f(x) x(i) = xi-h f2 = f(x) x(i) = xi g(i) = (f1-f2)/h2 end endfunction Habiendo cargado, mediante getf ..., el archivo que contiene estas dos funciones, su uso se puede hacer mediante ordenes semejantes a:

34 x = [3; 4]; gr = gradf(x)

CAP ITULO 4. PROGRAMAS

Si se desea una funcion que evalue el gradiente de varias funciones diferentes sin tener que editar cada vez la misma funcion f, entonces se debe pasar la funcion como uno de los parametros. //-------------------------------------------------------function fx = func1(x) fx = 3*x(1)^2 + 5*x(2) endfunction //-------------------------------------------------------function fx = func2(x) fx = 3*x(1) + 5*x(2)^3 endfunction //-------------------------------------------------------function g = grad(funcion, x) // gradiente de funcion en el punto x //******** h = 0.001 //******** h2 = 2*h n = max(size(x)) g = zeros(n,1) for i =1:n xi = x(i) x(i) = xi+h f1 = funcion(x) x(i) = xi-h f2 = funcion(x) x(i) = xi g(i) = (f1-f2)/h2 end endfunction //-------------------------------------------------------El llamado se puede hacer as : x = [3; 4], gr = grad(func1, x), fp = grad(func2, x)

4.6.2

Matriz escalonada reducida por las

En el siguiente ejemplo se obtiene la matriz escalonada reducida por las de una matriz dada. La unica diferencia con respecto al metodo tradicional que se ve en los cursos de

4.6. EJEMPLOS

35

Algebra Lineal es que se hace pivoteo parcial, es decir, se busca que el pivote sea lo mas grande posible en valor absoluto con el n de disminuir los errores de redondeo. Si la posicion del pivote es (i, j ), entonces se busca el mayor valor absoluto en las posiciones (i, j ), (i + 1, j ), ..., (m, j ) y se intercambia la la i con la la del mayor valor absoluto. Esta funcion merf permite escribir en pantalla resultados parciales de los calculos mediante la funcion escrMatr que a su vez llama la funcion escrVect. En merf se utiliza argn(2) que indica el numero de parametros en el llamado a la funcion y permite de nir los dos ultimos por defecto. El llamado merf(A) es equivalente a merf(A, 1.0e-12, 0) . El llamado merf(A, 1.0e-7) es equivalente a merf(A, 1.0e-7, 0). Para una matriz cualquiera, la orden norm(rref(a)-merf(a)) debe producir un valor muy cercano a cero. function escrVect(x, titulo) // escribe un vector en la pantalla

//********************* formato = %10.4f numNumerosPorLinea = 5 //********************* nnl = numNumerosPorLinea [n1, n2]= size(x) n = max(n1,n2) m = min(n1,n2) if m > 1 printf(ESCRVECT: x no es un vector.\n) end x = x(:) if argn(2) == 2 printf(%s :\n, titulo) end for i=1:n printf(formato, x(i)) if modulo(i, nnl) == 0 | i == n printf(\n) end end endfunction //------------------------------------------------------function escrMatr(A, titulo) // escribe una matriz en la pantalla

36

CAP ITULO 4. PROGRAMAS if argn(2) == 2 printf(%s :\n, titulo) end [m, n] = size(A) for i=1:m escrVect(A(i,:)) end endfunction //-------------------------------------------------function a = merf(A, eps0, IRP) // // // // // // Matriz escalonada reducida por filas. Se hace pivoteo parcial para buscar para disminuir los errores numericos. Si |t| <= eps0, se considera nulo. Si IRP = 1 se escriben resultados parciales. 0 no se escriben los resultados parciales

if argn(2) == 1 eps0 = 1.0e-12 IRP = 0 end if argn(2) == 2 IRP = 0 end if IRP == 1, escrMatr(a, matriz inicial), end [m, n] = size(a) eps0 = eps0*max(abs(A)) a = A i = 1 j = 1 while i <= m & j <= n [vmax, i0] = max(abs(a(i:m,j))) if vmax > eps0 imax = i+i0-1 if imax ~= i t = a(i,j:n) a(i,j:n) = a(imax,j:n) a(imax,j:n) = t

4.6. EJEMPLOS if IRP == 1 printf(interc. filas %d y %d\n, i, imax) escrMatr(a) end end aij0 = a(i,j) a(i,j+1:n) = a(i,j+1:n)/a(i,j) a(i,j) = 1 if IRP == 1 printf(dividir fila %d por el pivote: %f\n, i, aij0) escrMatr(a) end for k = 1:m if k ~= i a(k,j+1:n) = a(k,j+1:n) - a(k,j)*a(i,j+1:n) a(k,j) = 0.0 end end if IRP == 1 printf(buscar ceros en la columna %d\n, j) escrMatr(a) end i = i+1 else a(i:m,j) = zeros(m+1-i,1) end j = j+1 end endfunction

37

4.6.3

Aproximaci on polinomial por m nimos cuadrados

Dados m puntos (x1 , y1 ), (x2 , y2 ), ..., (xm , ym ), tales que xi = xj si i = j (no hay dos xi iguales), se desea encontrar un polinomio de grado n m 1 que pase lo mas cerca posible de estos puntos, en el sentido de m nimos cuadrados. Si n = m 1, el polinomio pasa exactamente por los m puntos. Los calculos que hay que hacer, utilizando las ecuaciones normales, son los siguientes: n 1 x 1 x2 x 1 1 1 x 2 x2 xn 2 2 A= , n 1 x m x2 m xm

38

CAP ITULO 4. PROGRAMAS

resolver (AT A) c = AT y,

donde,

c = c0 c1 c2 ... cn , p(x) = c0 + c1 x + c2 x2 + + cn xn , y ~ = Ac. Si n = m 1, no hay que resolver el sistema dado por las ecuaciones normales, se resuelve directamente, A c = y. function [p, yy, ind] = aproxPol(x, y, n) // Dados los vectore x, y, se desea // buscar p, el polinomio de aproximacion por // minimos cuadrados, tal que // p(x(i)) ~ y(i) // Si hay m puntos, se requiere que m >= n+1 // ind valdra 0 si hubo problemas // 1 si el proceso funciono bien. // yy contendra los valores del polinomio p en los // diferentes puntos x(i) p = poly([0], x, c) yy = [] ind = 0 x = x(:) y = y(:) [m, nx ] = size(x) [my, ny ] = size(y) if nx > 1 | ny > 1 printf(x, y deben ser vectores.\n) return end if m <> my

4.6. EJEMPLOS printf(x, y return end deben ser del mismo tama~no.\n)

39

if n < 0 | m < n+1 printf( n = %d inadecuado.\n, n) return end // no debe haber xi repetidos eps0 = 1.0e-16*max(abs(x)) for i=1:m-1 for j = i+1:m if abs(x(i)-x(j)) <= eps0 printf( x(i) repetidos.\n) return end end end ind = 1 A = zeros(m,n+1) for i=1:m A(i,1) = 1 xi = x(i) xij = 1 for j =1:n xij = xij*xi A(i,j+1) = xij end end if m == n+1 c = A\y else c = (A*A)\(A*y) end p = poly(c, x, c) yy = A*c endfunction

40

CAP ITULO 4. PROGRAMAS

4.6.4

Factores primos

En el siguiente ejemplo, para un entero n 2, se construye un vector la con los factores primos de n. En este vector la los factores primos aparecen repetidos segun la multiplicidad. Por ejemplo, para n = 72 = 23 32 , los factores seran 2, 2, 2, 3 y 3. La funcion factores hace uso de la funcion primerDiv que calcula el divisor mas peque~ no de un entero superior a uno. function [p, indic] = primerDiv(n) // calculo del primer divisor (divisor mas peque~no) // de un entero n >= 2. // indic valdra 0 si n es un dato malo, // 1 si el proceso se realizo bien. // Si n es primo, entonces p = n. p = 0 indic = 0 if n < 2 printf(PRIMERDIV: end if int(n) <> n printf(PRIMERDIV: end

n < 2.\n), return

no es entero.\n), return

ind = 1 n1 = floor(sqrt(n)) for p =2:n1 if modulo(n, p) == 0, return, end end p = n endfunction //-------------------------------------------------------function [p, indic] = factores(n) // Construccion del vector fila p con los factores // primos del entero n >= 2. // El producto de los elementos de p sera n. // indic valdra 0 si n es un dato malo, // 1 si el proceso se realizo bien. p = [] indic = 0

 4.7. MISCELANEA FACTORES:9e8.916.695A6 Tf0916.684AANEA 361.049e8.9924( )Tj -1. ..\n),10 0 9.123.1 0 0 cretur10 ANEAgravedad)10 MISCEL 5AMISCELMISCEL0 1 4628

42

CAP ITULO 4. PROGRAMAS // con velocidad inicial v0 en el tiempo if ~isdef(gravedad) gravedad = 9.8 end h = v0*t - 0.5*gravedad*t*t endfunction t.

Cuado se involucran vectores de manera adecuada, las comparaciones dan vectores de variables boolenas (valor %T o %F, correspondientes a 1 o 0). Por ejemplo, x = [ -2 -1.0e-10 y = [ 3 3 res = ( x <= y ) da como resultado el vector ! T T T F ! Las ordenes x = [ -2 -1.0e-10 -1.0e-20 res = ( abs(x) >= 1.0e-16 ) dan como resultado el vector ! T T F T ! La orden find permite obtener las posiciones de un vector booleano donde el valor es %T . Por ejemplo, las ordenes x = [ -2 -1.0e-10 -1.0e-20 4 ] p = find( abs(x) >= 1.0e-16 ) dan como resultado el vector ! 1. 2. 4. ! o sea, un vector con las posiciones de x donde hay valores no despreciables (ni nulos ni casi nulos). 4 ] -1.0e-20 3 4] 3]

Cap tulo 5 GRAFICAS

5.1 Dos dimensiones

Para hacer la gra ca de la funcion f : [a, b] R, basta con construir un vector con valores de x en el intervalo [a, b] y otro vector con los valores de f en los puntos del primer vector. Por ejemplo, a = -2; b = 3; x = a:0.01:b; y = sin(x); plot2d(x, y)

1.0 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1.0 2 1 0 1 2 3

hace la gra ca de seno(x) en el intervalo [2, 3]. 43

44

CAP ITULO 5. GRAFICAS

Obviamente los vectores x, y tienen el mismo tama~ no. Los valores de x se tomaron con un espaciamiento de 0.01. Esto hace que la gra ca se vea, no solo continua sino tambien suave. Si el espaciamiento es muy grande, por ejemplo, a = -2; b = 3; x = a:0.5:b; y = sin(x); plot2d(x, y) dara una gra ca continua pero poligonal y no suave. De hecho siempre Scilab hace l neas poligonales, pero si el espaciamiento es muy peque~ no la l nea poligonal se confunde con una curva suave. En resumen, el espaciamiento debe ser peque~ no, pero es ino cioso hacerlo muy peque~ no. Las ordenes siguientes permiten obtener la gra ca de una parabola: a = -2; b = 3; x = a:0.05:b; y = 0.5*x.*x - 2; plot2d(x, y) En una sesion de Scilab, la primera vez que se hace una gra ca, esta aparece inmediatamente en la pantalla. Cuando se da la orden para una segunda gra ca, esta es creada pero no aparece automaticamente en la pantalla. Es necesario, mediante un clic, activar la ventana de la gra ca. Muy posiblemente despues de la segunda orden plot2d, en la gra ca apareceran las dos \curvas" superpuestas. Para limpiar la ventana gra ca se utiliza xbasc(). Observe la diferencia entre los dos grupos de ordenes siguientes: plot2d(x,sin(x)) plot2d(x,cos(x))

plot2d(x,sin(x)) xbasc(), plot2d(x,cos(x)) Si se da la orden plot(x, y, t, sen(t), Ejemplo 1); en la gra ca apareceran ademas 3 etiquetas (o rotulos), la primera en el eje x, la segunda en el eje y y la tercera sera la etiqueta general de la gra ca. En la misma gura pueden aparecer varias funciones. Para esto, los datos de plot2d deben estar en columnas o en matrices. Si x, y, z son vectores columna con el mismo numero de las, entonces se puede dar la orden

5.1. DOS DIMENSIONES plot2d(x, [y z] )

45

cuyo efecto es tener en la misma gura las gra cas producidas por las ordenes plot2d(x,y) y plot2d(x,z) . Por ejemplo, x = (-2:0.05:3); y = sin(x); z = cos(x); plot2d(x, [y z sin(2*x)]) En este caso Scilab asigna colores diferentes a cada una de las curvas. Para conocer la correspondencia entre colores y curvas se pueden colocar letreros (\legend") para las curvas. Se utiliza la opcion leg separando los letreros por medio de @. plot2d(x, [y z sin(2*x)], leg="sen(x)@cos(x)@sen(2x)" ) Despues de haber hecho una gra ca, se le pueden agregar letreros. Se usa xtitle , que tiene 3 parametros, todos los tres deben ser cadenas. El primero para el letrero general, el segundo para el eje horizontal y el tercero para el eje vertical. Por ejemplo xtitle("Funciones trigonometricas", "x", "") En resumen, una de las formas de utilizar plot2d es plot2d(x, a), donde x es un vector columna y a es una matriz con el mismo numero de las que x. En la gura habra tantas gra cas como columnas tenga la matriz a. Una forma mas general es plot2d(u, v), donde u y v son matrices del mismo tama~ no, de m las y n columnas. Entonces en la misma gura apareceran la gra ca de la primera columna de u con la primera columna de v, la gra ca de la segunda columna de u con la segunda columna de v, etc. En la gura producida por el siguiente ejemplo esta la gra ca de seno entre 0 y 3.14 y la gra ca de la tangente entre 1 y 1. n = 100; a = 0; b = 3.14; h = (b-a)/n; x1 = (a:h:b); y = sin(x1); a = -1; b = 1; h = (b-a)/n; x2 = (a:h:b); z = tan(x2); plot2d([x1 x2], [y z]); Tambien se puede, usando fplot2d , obtener la gra ca de una funcion f de nida mediante una funcion de Scilab. En este caso, solo se necesita el vector de valores de xi . Este vector debe ser monotono, es decir, creciente o decreciente. Supongamos que se ha de nido la siguiente funcion function fx = func4(x) fx = sin(x)-tan(x)

46 endfunction

CAP ITULO 5. GRAFICAS

y que se ha cargado el archivo donde esta de nida, mediante getf ... Entonces se puede obtener la gra ca mediante u = (-1:0.01:1); fplot2d(u, func4) Si la funcion es sencilla, no es necesario crear la funcion en un archivo. Se puede obtener la gra ca mediante: t = (-1:0.01:1.2); deff("[y] = func4b(x)", "y = x.*x + 2.5") fplot2d(t, func4b) Algunas veces se desea, no obtener la curva, sino dibujar los puntos, sin que haya \continuidad", o sea, que los puntos no esten unidos por un segmento de recta. Una forma de hacerlo es mediante un argumento adicional, el estilo. Las siguientes ordenes muestran un ejemplo. xbasc(); x = (-3.2:0.2:3.2); y = sin(x); tipo = -1; plot2d( x, y, tipo) Si el estilo, en este ejemplo la variable tipo, es menor o igual a cero, en la gra ca aparecen los puntos representados por un s mbolo. Los valores pueden ser 0, 1, 2, ..., 9 El signi cado de los cinco primeros es el siguiente: 0: -1 : -2 : -3 : -4 : Obviamente en la gra ca aparecen reducciones de estos s mbolos. Es posible dibujar, al mismo tiempo, dos listas de puntos (con igual cantidad de puntos), cada lista con un s mbolo diferente. Por ejemplo: xbasc(); x = (-3.2:0.2:3.2); plot2d( [x x], [sin(x) cos(x)], [-1 -2] ) +

5.2. TRES DIMENSIONES

47

5.2

Tres dimensiones

Sea f : R2 R, por ejemplo, f (x, y ) = 5x2 y 2 . Para hacer la gra ca de la funcion restringida al rectangulo [a, b] [c, d], basta con construir un vector con valores de la primera variable en el intervalo [a, b], otro con valores de la segunda variable en el intervalo [c, d] y una matriz con los valores f (..., ...). Por ejemplo, u = v = m = w = for (-2:0.05:2); (-3:0.1:3); size(u,1); n = size(v,1); zeros(m,n); i=1:m for j = 1:n w(i,j) = 5*u(i)^2 - v(j)^2; end

end plot3d(u, v, w); De manera analoga a dos dimensiones, es posible gra car por medio de fplot3d una funcion f de nida por medio de una funcion de Scilab. Basta con dar dos vectores, uno con las coordenadas de la primera variable y otro con las coordenadas de la segunda variable. Supongamos que se ha de nido la siguiente funcion function fst = func5(s, t) fst = 5*s^2 - t^2 endfunction y que se ha cargado el archivo donde esta de nida, mediante getf ... Entonces se puede obtener la gra ca mediante u = (-4:0.05:4); v = (-5:0.05:3); fplot3d(u, v, func5) Esta funcion fplot3d construye la matriz de valores f (u[i],v[j]), de tama~ no, para este ejemplo, 161 221. Esto se demora algunos segundos. Ademas si u y v son muy grandes, la matriz podr a ser demasiado grande. En este caso Scilab muestra un aviso semejante a: stack size exceeded (Use stacksize function to increase it)

48

CAP ITULO 5. GRAFICAS

indicando que el tama~ no de la pila no es su ciente para lo solicitado. Este tama~ no se puede modi car mediante la funcion stacksize, por ejemplo, stacksize(2000000) Se pueden obtener curvas de nivel mediante la funcion contour. Su uso es muy semejante al de plot3d pero con un parametro adicional que indica el numero de curvas de nivel. Por ejemplo: u = v = m = w = for (-2:0.05:2); (-3:0.1:3); size(u,1); n = size(v,1); zeros(m,n); i=1:m for j = 1:n w(i,j) = 5*u(i)^2 + v(j)^2; end

end contour(u, v, w, 5); Tambien, mediante fcontour, se pueden obtener curvas de nivel de una funcion de nida en una funcion de Scilab. Su uso es analogo al de fplot3d pero con un parametro adicional que indica el numero de curvas de nivel. Por ejemplo: u = (-4:0.2:4); v = (-4:0.2:5); fcontour(u, v, func5, 7)

5.3

Otras funciones

La siguiente lista contiene los nombres de la mayor a de la funciones para gra cas. addcolor alufunctions black bode champ champ1 chart colormap contour contour2d contour2di contourf dragrect drawaxis driver edit_curv errbar eval3d eval3dp evans fac3d fchamp fcontour fcontour2d fec fgrayplot fplot2d fplot3d fplot3d1 gainplot genfac3d geom3d getcolor getfont getlinestyle getmark getsymbol gr_menu graduate graycolormap grayplot graypolarplot hist3d histplot hotcolormap isoview legends locate m_circle Matplot Matplot1 milk_drop nf3d nyquist param3d param3d1 paramfplot2d plot plot2d plot2d1 plot2d2 plot2d3 plot2d4 plot3d plot3d1 plot3d2 plot3d3 plotframe plzr polarplot printing

 DE UN ARCHIVO POSTSCRIPT 5.4. CREACION replot rotate scaling sd2sci secto3d Sfgrayplot Sgrayplot sgrid square subplot titlepage winsid xarc xarcs xarrows xaxis xbasc xbasimp xbasr xchange xclea xclear xclick xclip xdel xend xfarc xfarcs xfpoly xfpolys xfrect xget xgetech xgetmouse xgraduate xgrid xinfo xinit xlfont xload xname xnumb xpause xpoly xpolys xrect xrects xrpoly xs2fig xsave xsegs xselect xset xsetech xsetm xstring xstringb xstringl xtape xtitle zgrid

49

5.4

Creaci on de un archivo Postscript

Una vez hecha la gra ca, la creacion de archivos Postscript o de otros tipos de formato puede ser realizada por medio de la barra de menu de la ventana de la gra ca, con las opciones File Export Postscript nombre-del-archivo El nombre del archivo debe ser dado sin extension. Observe que para exportar hay otras posibilidades: Postscript No Preamble: Postscript-Latex Xfig GIF PPM Otra forma de crear archivos postscript es mediante las ordenes: driver xinit xend. Utilice la ayuda de Scilab par obtener informacion sobre su uso. La orden xbasimp permite enviar una gra ca a una impresora postscript o a un archivo.

Cap tulo 6 METODOS NUMERICOS

6.1 Sistemas de ecuaciones lineales

Como se vio anteriormente, si se tiene un matriz A de tama~ no n n y un vector b de tama~ no n 1, la manera e ciente de reolver el sistema Ax = b, es mediante x = A\b Si la matriz es singular o casi singular, Scilab muestra una advertencia. Si la matriz es de nida positiva, resolver el sistema por medio de x = A\b funciona bien, pero es mas e ciente resolverlo usando la factorizacion de Cholesky: y = zeros(n,1); x = zeros(n,1); u = chol(a); y(1) = b(1)/u(1,1); for i=2:n y(i) = ( b(i) - u(1:i-1,i)*y(1:i-1) )/u(i,i); end x(n) = y(n)/u(n,n); for i=n-1:-1:1 x(i) = ( y(i) - u(i,i+1:n)*x(i+1:n) )/u(i,i); end La diferencia en el tiempo se nota para valores grandes de n. Por ejemplo, si n = 1500, la solucion x = A\b gasta 12.33 \segundos" (en realidadson unidades de tiempo, no exactamente segundos, en un computador espec co). En cambio, obtener la solucion utilizando la factorizacion de Cholesky dura 2.78 segundos. Para medir el tiempo se puede puede utilizar la funcion timer(). 50

 POR M 6.2. SOLUCION INIMOS CUADRADOS

51

Para economizar memoria se hubiera podido efectuar el calculo de y y de x sobreescribiendo encima de b. Una version, con menos ordenes, pero ine ciente, podr a ser: u = chol(a); y = (u\b); x = u\y; Por ejemplo, con n = 1500, estas tres ordenes duran 15.12 segundos. Una matriz dispersa (sparse ) es aquella que tiene muchos ceros y muy pocos elementos no nulos. Si 8 2 A= 0 0 0 0 0 0 1.12 5 0.5 0 0 0 10 3.14 0 0 0 4 1 0 0 0 20

basta con almacenar, para cada entrada no nula, el ndice de la la el ndice de la columna y el valor. Entonces en Scilab se requiere una matriz de dos columnas con los ndices (numeros enteros entre 1 y n) y un vector con los valores: ij = [1 1; 1 5; 2 1 ; 2 2; 2 3;3 3; 3 4; 4 4; 4 5; 5 5] v = [-8 1.12 -2 5 0.5 10 3.14 4 -1 20] A = sparse(ij, v) b = [-6.88 3.5 13.14 3 20] x = A\b Construida de esa manera, se puede resolver el sistema Ax = b utilizando A\b. Esta forma de almacenamiento es util cuando A es muy grande y muy dispersa.

6.2

Soluci on por m nimos cuadrados

Cuando se tiene el sistema Ax = b, donde A es un matriz m n, m n, es muy probable que el sistema no tenga solucion, dado que hay mas ecuaciones que incognitas. En este caso se busca obtener un x que minimice el cuadrado de la norma euclidiana del error, o sea

52

CAP ITULO 6. METODOS NUMERICOS

minimizar ||Ax b||2 2 Por ejemplo, A = [1 2; 3 4; 5 6] b = [3 7 11.1] En este caso, la solucion por m nimos cuadrados se puede obtener mediante x = A\b o mediante x = pinv(A)*b La primera forma es mas e ciente.

6.3

Una ecuaci on no lineal

Para resolver f (x) = 0, donde f es una funcion de variable y valor real, se utiliza fsolve. Por ejemplo para resolver x ex cos(x) + 0.1 = 0, 1 + x2 es necesario de nir una funcion de Scilab donde este f y despues utilizar fsolve. function fx = func156(x) fx = ( x - exp(x) )/( 1 + x*x ) - cos(x) + 0.1 endfunction Despues de haber cargado esta funcion, se utiliza fsolve dandole como parametros, la aproximacion inicial y la funcion: r = fsolve(0, func156) Con otra aproximacion inicial podr a dar otra ra z. Un parametro opcional, que puede acelerar la obtencion de la solucion, es otra funcion de Scilab donde este de nida la derivada.

6.4. SISTEMA DE ECUACIONES NO LINEALES

53

6.4

Sistema de ecuaciones no lineales

Dada una funcion f : Rn Rn se desea encontrar un vector x tal que f (x) = 0. Aqu 0 es el vector nulo. Por ejemplo, x2 1 + x 1 x2 + x 3 3 = 0 , 2x1 + 3x2 x3 5 = 0 , (x1 + x2 + x3 )2 10x3 + 1 = 0 . La funcion f debe estar de nida en una funcion de Scilab, donde, dado un vector columna x, calcule el vector columna f (x). Una vez escrita la funcion y cargado el archivo donde esta f , se utiliza fsolve . Para el ejemplo anterior la funcion puede ser function fx = f3(x) n = size(x,1) fx = zeros(n,1) fx(1) = x(1)^2 + x(1)*x(2) + x(3) - 3 fx(2) = 2*x(1) + 3*x(2)*x(3) - 5 fx(3) = (x(1) + x(2) + x(3))^2 - 10*x(3) + 1 endfunction El llamado de fsolve es de la forma [x, fx, info] = fsolve([3 4 5], f3) El primer parametro de entrada es la aproximacion inicial de la solucion. El segundo es el nombre de la funcion donde esta de nida f . Hay tres parametros de salida. Unicamente el primero es obligatorio. Es un vector con la solucion, si esta fue obtenida. El vector fx es el vector f (x). El parametro info puede tener 5 valores: 0 : numero inadecuado de parametros de entrada. 1 : se obtuvo la solucion con la precision deseada. 2 : muchos llamados a la funcion. 3 : la tolerancia permitida es muy peque~ na. 4 : el metodo avanza muy lentamente. Si se dispone de una funcion donde esta de nida la matriz jacobiana, entonces fsolve trabajara mucho mejor. Por ejemplo, function J = jacf3(x) t = 2*(x(1)+x(2)+x(3)) J = [ 2*x(1)+x(2) x(1)

54 2 t endfunction 3*x(3) t

CAP ITULO 6. METODOS NUMERICOS 3*x(2) t-10 ]

Entonces el llamado de fsolve es de la forma [x, fx, info] = fsolve([3 4 5], f3, jacf3) El problema de encontrar la solucion, por m nimos cuadrados, de un sistema de ecuaciones no lineales esta en el cap tulo siguiente.

6.5

Integraci on num erica

Para obtener una aproximacion del valor de una integral de nida, por ejemplo,
0.5 0.1

ex dx

se utiliza intg . Para eso es necesario de nir en Scilab la funcion que se va a integrar. Puede ser, directamente en el ambiente Scilab: deff([y] = f53(x), y = exp(-x*x)) I = intg(0.1, 0.5, f53) Tambien se puede de nir una funcion en un archivo .sci function fx = f57(x) fx = exp(-x*x) endfunction y despues de cargarla, dar la orden I = intg(0.1, 0.5, f57) Algunas veces no se conoce una expresion de la funcion f , pero se conoce una tabla de valores (xi , f (xi )), o simplemente una tabla de valores (xi , yi ). Supongamos, as lo requiere Scilab, que la lista de valores (x1 , y1 ), ..., (xn , yn ) esta ordenada de manera creciente de acuerdo a los xi , o sea, x1 < x2 < < xn .

Para obtener el valor aproximado de la integral, entre x1 y xn , de la funcion f (representada por los valores (xi , yi )), es necesario tener dos vectores con los valor xi y yi , y utilizar la funcion inttrap . x = [0.1 0.15 0.2 0.25 0.3 0.4 0.5] y = [ 0.9900 0.9778 0.9608 0.9394 0.9139 0.8521 0.7788] I = inttrap(x, y)

 NUMERICA 6.6. DERIVACION

55

6.6

Derivaci on num erica

Para obtener una aproximacion numerica de la derivada de una funcion en un punto se puede usar derivative . Si f57 es la funcion de nida anteriormente, una vez cargada, la orden d = derivative(f57, 0.3) da una aproximacion numerica de la derivada de la funcion en x = 0.3. Si se tiene una funcion f : Rn R, por ejemplo, function fx = func90(x) // funcion de Rosenbrock t1 = x(2)-x(1)*x(1) t2 = 1-x(1) fx = 10*t1*t1 + t2*t2 + 4 endfunction entonces la orden derivative(func90, [3 4]) da como resultado un vector la con el gradiente de f en x = (3, 4). El segundo parametro, el punto x, debe ser un vector columna. Si se tiene una funcion f : Rn Rm , entonces f (x) = ( f1 (x), f2 (x), ..., fm (x) ), donde fi : Rn R . Por ejemplo, function Fx = func157(x) Fx = zeros(4,1); Fx(1) = x(1)^2 + x(1)*x(2) + x(3) - 3; Fx(2) = 2*x(1) + 3*x(2)*x(3) - 5; Fx(3) = ( x(1) + x(2) + x(3) )^2 - 10*x(3) + 1; Fx(4) = 10*( x(1) + x(2) + x(3) - 3.3 ); endfunction entonces la orden derivative(func157, [3 4 5]) da como resultado la matriz jacobiana de la funcion f (de nida en f157 ) evaluada en el punto x = (3, 4, 5). El segundo parametro, el punto x, debe ser un vector columna. La matriz jacobiana es de tama~ no m n (en este ejemplo: 4 3 ). En la la i esta el gradiente de fi .

56

CAP ITULO 6. METODOS NUMERICOS

6.7

Interpolaci on

Cuando hay m puntos (x1 , y1 ), (x2 , y2 ), ..., (xm , ym ) se desea obtener la funcion interpolante, una funcion que pase por esos puntos, con el objetivo de evaluarla en otros valores x intermedios. La funcion interpln permite hacer interpolacion lineal. Tiene dos parametros, el primero es una matriz de dos las. La primera la tiene los valores xi . Deben estar en orden creciente. La segunda la tiene los valores yi . El segundo parametro es un vector donde estan los valores x en los que se desea evaluar la funcion interpolante (af n por trozos). x = [.4trozos p8 p8 rozos). rozos). (a Td )Tj 14.0395 0 Td (p6trozos).)Tj /R27 117

 POR M 6.8. APROXIMACION INIMOS CUADRADOS

57

una matriz de tama~ no m (n + 1), resolver por m nimos cuadrados un sistema de ecuaciones. Generalmente, despues se evalua este polinomio en algunos valores (pueden ser los mismos xi u otros). x y m t n = = = = = [0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2]; [4.8 6.2 6.8 7.2 7.8 9.2 8.8 9.2 8.8 9.2 7.8]; size(x,1); (x(1):0.01:x(m)); 2;

A = zeros(m, n+1); for i=0:n A(:,i+1) = x.^i; end cf = A\y; p = poly(cf, x, c); ft = horner(p, t); xbasc() plot2d(t, ft) Para la aproximacion no es indispensable que los xi esten ordenados de manera creciente, pero en este ejemplo es necesario para la gra ca. Tambien se puede hacer, de manera analoga, aproximacion por m nimos cuadrados utilizando una combinacion lineal de otras funciones (trigonometricas, exponenciales, logar tmicas...). Cambia simplemente la manera de construir la matriz A y la manera de evaluar la funcion en los valores del vector t. Scilab tiene la funcion lsq_splin que permite aproximacion por m nimos cuadrados con funciones spline. Los datos son: los valores xi , ordenados de manera creciente, los valores yi , y los puntos de ruptura (breakpoint ) ordenados de manera creciente. Con ellos se determina la funcion spline. Despues se evalua en los valores deseados. x = [0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2]; y = [4.8 6.2 6.8 7.2 7.8 9.2 8.8 9.2 8.8 9.2 7.8]; nx = size(x,1); t = (x(1):0.01:x(nx)); X = [0 1 2] [Y, d] = lsq_splin(x, y, X); ft = interp(t, X, Y, d); xbasc() plot2d(t, ft)

58

CAP ITULO 6. METODOS NUMERICOS

6.9

Ecuaciones diferenciales ordinarias

Una ecuacion diferencial de primer orden con condicion inicial se escribe, generalmente, en la forma

y = f (x, y ), y (x0 ) = y0 . Por ejemplo, x+y +4, x2 + y 2 y (2) = 3. y = Antes de utilizar la funcion ode , es necesario crear en Scilab la funcion f y cargarla. La funcion ode evalua aproximaciones del valor de y en valores t1 , t2 , ..., tp . Despues de de nir y cargar function fxy = func158(x, y) fxy = ( x + y )/( x*x + y*y ) + 4 endfunction se obtiene la solucion aproximada mediante x0 = 2 y0 = 3 t = [2 2.1 2.2 2.5 2.9] yt = ode(y0, x0, t, func158)

6.10

Sistema de ecuaciones diferenciales

Un sistema de ecuaciones diferenciales ordinarias de primer orden con condiciones iniciales se escribe, generalmente, en la forma

y = f (x, y ), y (x0 ) = y0 .

6.10. SISTEMA DE ECUACIONES DIFERENCIALES pero ahora, y , y y y0 son vectores en Rn ; f : R1+n Rn . Por ejemplo, 2y1 + x3 y2 , x 3 y2 = y2 x y1 (1) = 1 y2 (1) = 1. y1 = Despues de de nir y cargar function fxy = func43(x, y) fxy = zeros(2,1) fxy(1) = 2*y(1)/x + x^3*y(2) fxy(2) = -3*y(2)/x endfunction se utiliza la misma funcion ode , pero con los parametros de dimension adecuada. x0 = 1 y0 = [-1 1] t = (1:0.2:2) yt = ode(y0, x0, t, func43)

59

En este caso, yt es un matriz de dos las. En la la i estan las aproximaciones de los valores de yi (tj ).

Cap tulo 7 OPTIMIZACION

7.1 Optimizaci on lineal

Hay varias maneras o formas de presentacion de problemas de OL, optimizacion lineal (o programacion lineal). Entre las de uso mas frecuente estan la forma canonica y la forma estandar. La forma que mas rapidamente se adapta a la solucion en Scilab, mediante la funcion linpro , es la siguiente: min z = cT x Ai x = bi , i = 1, ..., m0 Ai x bi , i = m0 + 1, ..., m u x v.

(7.1)

Obviamente, siempre se puede pasar de una forma a otra sin di cultad. Para el caso de Scilab, todos los problemas deben ser de minimizacion y las desigualdades, si las hay, diferentes de las restricciones de caja (u x v ), deben ser de la forma ai1 x1 + ai1 x2 + + ai1 xn bi . La matriz A es de tama~ no m n; los vectores c, u, v Rn1 , b Rm1 (todos son vectores columna). Las igualdades, si las hay, deben estar antes de las desigualdades. Frecuentemente, en los problemas de OL, todas las variables deben ser no negativas. En Scilab esto se debe considerar como restricciones de caja o como restricciones usuales.

7.1.1

Desigualdades

El problema mas sencillo para resolver por medio de linpro es el siguiente: min cT x 60

 LINEAL 7.1. OPTIMIZACION Ax b. La orden en Scilab es: [x, lagr, z] = linpro(c, A, b)

61

El signi cado de los tres resultados es el siguiente : x es el vector solucion; lagr es un vector con los coe cientes de Lagrange; z es el valor optimo de la funcion objetivo, es decir, cT x. Para resolver el siguiente problema min x1 1.4x2 x1 + x2 400 x1 + 2x2 580 x1 300 x 0, es necesario reescribir las restricciones de no negatividad: min x1 1.4x2 x1 + x 2 x1 + 2x2 x1 x1 x2 Su solucion mediante Scilab: c = A = b = [x, [-1 -1.4 ]; [ 1 1; 1 2; 1 0; -1 0 ; 0 -1]; [400 580 300 0 0]; lg, z ] = linpro(c, A, b) 400 580 300 0 0.

(7.2)

7.1.2

Desigualdades y restricciones de caja

Se usa la orden [x, lg, z ] = linpro(c, A, b, u, v) para resolver min cT x

62 Ax b u x v.

CAP ITULO 7. OPTIMIZACION

El problema (7.2) se puede resolver con las siguientes ordenes: c = A = b = [x, [-1 -1.4 ]; [ 1 1; 1 2; 1 0]; [400 580 300]; lg, z ] = linpro(c, A, b, zeros(2,1), [])

Es interesante observar que como no hay restricciones superiores de caja, entonces se utiliza un vector \vac o". En el problema (7.2) la restriccion x1 300 se puede considerar como una restriccion de caja.

min x1 1.4x2 x1 + x2 400 x1 + 2x2 580 x1 300 x 0.

As planteado, se puede resolver mediante las siguientes ordenes: c = A = b = v = [x, [-1 -1.4 ]; [ 1 1; 1 2]; [400 580]; [ 300 1.0e90]; lg, z ] = linpro(c, A, b, zeros(2,1), v)

7.1.3

Igualdades, desigualdades y restricciones de caja

Es el caso general (7.1) planteado al comienzo del cap tulo. Este problema se resuelve mediante la orden [x, lg, z ] = linpro(c, A, b, u, v, m0) El siguiente problema min x1 1.4x2 x1 + x2 = 400

 CUADRATICA 7.2. OPTIMIZACION x1 + 2x2 580 x1 300 x 0,

63

se resuelve mediante: c = A = b = v = [x, [-1 -1.4 ]; [ 1 1; 1 2]; [400 580]; [ 300 1.0e90]; lg, z ] = linpro(c, A, b, zeros(2,1), v, 1)

7.2

Optimizaci on cuadr atica

El problema general de optimizacion cuadratica es muy parecido al caso general de OL. La unica diferencia consiste en que la funcion objetivo tiene una parte cuadratica 1 T x Qx. Las restricciones siguen exactamente el mismo esquema. 2 1 min q (x) = xT Hx+cT x 2 Ai x = bi , i = 1, ..., m0 Ai x bi , i = m0 + 1, ..., m u x v. La matriz H debe ser simetrica. Se usa la funcion quapro . Su llamado y uso son muy semejantes a los de la funcion linpro. [x, lg, qx ] = quapro(H, c, A, b, u, v, m0) Los casos de desigualdades y, desigualdades y restricciones de caja, son completamente analogos a los de OL. Consideremos el siguiente problema:

min (x1 3)2 + (x2 4)2 x1 + 2x2 3 3x1 + x2 4 x 0.

64 Entonces: H=

CAP ITULO 7. OPTIMIZACION 6 2 0 . , c= 8 0 2

Su solucion se puede obtener mediante las ordenes H = c = A = b = [x,

[ 2 0; 0 2]; [-6; -8]; [ 1 2; 3 1]; [ 3 ; 4]; lagr, qx ] = quapro(H, c, A, b, zeros(2,1), [])

7.3
7.3.1

Optimizaci on no lineal
Optimizaci on no restringida

El problema que se trata de resolver es simplemente minimizar una funcion f : Rn R, min f (x) Para su solucion en Scilab se usa la funcion optim. Se necesita haber de nido previamente una funcion en Scilab donde se calcule f (x) y el gradiente f (x). Esta funcion que por el momento llamaremos f_obj debe tener la siguiente forma [fx, gr, ind2] = f_obj(x, ind1) Parametros de entrada: x: punto donde se va a evaluar f (x); si ind1 = 2, se evalua f (x); si ind1 = 3, se evalua f (x); si ind1 = 4, se evaluan f (x) y f (x). Parametros de salida: fx: valor f (x); gr: gradiente f (x); si ind2 < 0, no se puede evaluar f (x); si ind2 = 0, se interrumpe la optimizacion. Consideremos dos problemas, el primero esta dado por la llamada funcion de Rosenbrock. Las soluciones de estos problemas son x = (1, 1) y x = (3, 4).
2 2 min f (x1 , x2 ) = 10(x2 x2 1 ) + (1 x1 ) , min f (x1 , x2 ) = (x1 3)2 + (x2 4)2 .

Antes de resolver estos dos problemas de minimizacion mediante optim es necesario construir un archivo, por ejemplo de nombre func.sci

 NO LINEAL 7.3. OPTIMIZACION function [fx, gr, ind2] = rosenbr(x, ind) // funcion de Rosenbrock // f(x) = 10(x2-x1^2) + (1-x1)^2 fx = 0.0; gr = 0.0*x; ind2 = 1; if ind == 2 | ind == 4 fx = 10*(x(2)-x(1)^2)^2 + (1-x(1))^2; end if ind == 3 | ind == 4 gr(1,1) = -40*x(1)*(x(2)-x(1)^2) - 2*(1-x(1)); gr(2,1) = 20*(x(2)-x(1)^2); end endfunction //-----------------------------------------------function [fx, gr, ind2] = ensayo(x, ind) fx = 0.0; gr = 0.0*x; ind2 = 1; if ind == 2 | ind == 4 fx = (x(1)-3)^2 + (x(2)-4)^2; end if ind == 3 | ind == 4 gr(1,1) = 2*x(1)-6; gr(2,1) = 2*x(2)-8; end endfunction La funcion optim se puede llamar de la siguiente forma: [fx, x] = optim(f_obj, x0) Para nuestros dos ejemplos espec cos, las ordenes pueden ser getf c:\...\func.sci [fxopt, xopt] = optim(rosenbr, [5; 6]) [fxopt, xopt] = optim(ensayo, [0; 0])

65

La funcion optim admite escoger el metodo de minimizacion entre cuasi Newton, gradiente conjugado y no diferenciable. Para ello es necesario utilizar qn gc o nd. La opcion por defecto es qn. Por ejemplo, [fx, x]=optim(rosenbr, x0, qn) [fx, x]=optim(rosenbr, x0, gc) [fx, x]=optim(rosenbr, x0, nd)

66

CAP ITULO 7. OPTIMIZACION

7.3.2

Restricciones de caja

La funcion optim tambien admite restricciones de caja, o sea, sirve para un problema de la forma
2 2 min f (x1 , x2 ) = 10(x2 x2 1 ) + (1 x1 ) 3 x1 5 2 x2 10 .

La ordenes pueden ser: u = [3; 2]; v = [5; 10]; x0 = [0; 0]; [fx, x] = optim(rosenbr, b, u, v, x0) [fx, x] = optim(rosenbr, b, u, [fx, x] = optim(rosenbr, b, u, v, x0, qn) v, x0, gc)

El parametro b (bound ) indica que los dos vectores columna siguientes son las cotas de caja. Tambien es posible escoger el criterio de parada, para ello se usa el parametro ar . Ver la ayuda de la funcion optim.

7.4

M nimos cuadrados no lineales

Dada una funcion f : Rn Rm , (usualmente n m), se desea minimizar el cuadrado de la norma euclidiana de f (x), min ||f (x)||2 . Por ejemplo,
2 f (x1 , x2 , x3 ) = (x2 1 + x1 x2 + x3 3, 2x1 + 3x2 x3 5, (x1 + x2 + x3 ) 10x3 + 1, 10(x1 + x2 + x3 ) 3.3).

La solucion se hace mediante leastsq . De manera analoga al caso de sistemas de ecuaciones no lineales, es necesario escribir una funcion de Scilab donde se de ne f (dado un vector columna x calcula el vector columna f (x) ). Por ejemplo function fx = fmc(x) m = 4;

7.4. M INIMOS CUADRADOS NO LINEALES fx = zeros(m,1); fx(1) = x(1)^2+x(1)*x(2)+x(3)-3; fx(2) = 2*x(1)+3*x(2)*x(3)-5; fx(3) = (x(1)+x(2)+x(3))^2-10*x(3)+1; fx(4) = 10*(x(1)+x(2)+x(3)-3.3); endfunction El llamado de leastsq es muy parecido al de fsolve. Por ejemplo, [fx, x] = leastsq(fmc, [3 5 7])

67

Tambien el metodo trabaja mejor si se tiene una funcion de Scilab donde se calcula la matriz con las derivadas parciales. Por ejemplo, function D = derfmc(x) t = 2*(x(1)+x(2)+x(3)); D = [ 2*x(1)+x(2) x(1) 2 3*x(3) t t 10 10 endfunction La orden pueden ser: [fx, x] = leastsq(fmc, derfmc, [4 5 6])

1 3*x(2) t-10 10 ];

7.4.1

M nimos cuadrados con restricciones de caja

En Scilab se puede resolver el problema de m nimos cuadrados no lineales con restricciones de caja: min ||f (x)||2 u x v. Las ordenes pueden ser: u = [1.2 0 1.1]; v = [10 10 10]; [fx, x] = leastsq(fmc, b, u, v, [4 5 6]) o para mayor e ciencia, [fx, x] = leastsq(fmc, derfmc, b, u, v, [4 5 6])

68 Otras herramientas de optimizacion son: lmisolver, lmitool: linear matrix inequalities, karmarkar : metodo de Karmarkar. para OL. semidef : optimizacion semide nida,

CAP ITULO 7. OPTIMIZACION

 INDICE ANAL ITICO

, 9 b, 66 c, 10 coeff, 10 e, 9 gc, 65 nd, 65 qn, 65 r, 10 roots, 10 v, 9 (:), 15 .sce, 22 .sci, 24 \, 20 %, 9 timer, 50 abort, 31 abs, 6, 9 acos, 6 acosh, 6 aleatorio, 7 algoritmo de Euclides, 31 anidamiento, 29 ans, 5 apropos, 8 aproximacion polinomial por m nimos cuadrados, 37 aproximacion por m nimos cuadrados, 56 archivo Postscript, 49 arcocoseno, 6 arcocoseno hiperbolico, 6 arcoseno, 6 arcoseno hiperbolico, 6 arcotangente, 6 69 arcotangente hiperbolica, 7 asignacion, 3 asin, 6 asinh, 6 atan, 6 atanh, 7 ayuda, 8 b, 66 backslash, 20 barra de menu, 8, 23 break, 31 C, 32 c, 10 C++, 24 carpeta actual, 26 por defecto, 26 case, 30 ceil, 7, 8 chdir, 27 chol, 19 Cholesky, 19, 50 coeff, 10 coe ciente, 10 columna, 13, 14 comentarios, 24 comparaciones, 27 complejos, 9 conj, 9 conjugado, 9 constantes prede nidas, 9 contour, 48 control de ujo, 27 cos, 7

70 coseno, 7 coseno hiperbolico, 7 cosh, 7 cotangente, 7 cotangente hiperbolica, 7 cotg, 7 coth, 7 cuasi Newton,, 65 curvas de nivel, 48 def, 54 derivacion numerica, 55 derivadas parciales, 67 derivative, 55 descomposicion en valores singulares, 20 desigualdades, 60{62 desigualdades lineales de matrices, 68 desviacion estandar, 18 det, 16 determinante, 16 diag, 17 diagonal de una matriz, 17 Dialog, 8 diferente, 27 disp, 32 dispersa, 51 division, 6, 9 elemento a elemento, 15 e, 9 ecuacion no lineal, 52 ecuaciones diferenciales ordinarias, 58 ecuaciones no lineales, 53 else, 28, 30 end, 28{30 endfunction, 24 ENPC, 1 enteros, 4 espacios en blanco, 5 espectro, 19 estructuras de control, 26, 27 evaluar un polinomio, 11 Exec, 23

INDICE ANAL ITICO exec, 23, 24, 27 exists, 41 exp, 7 expm, 17 exponencial de una matriz, 17 Export, 49 eye, 14 factores primos, 40 factorizacion de Cholesky, 19 LU, 20 QR, 20 fcontour, 48 la, 14 File, 23, 25, 49 find, 42 fix, 7 floor, 7, 8 for, 29 format, 9 formato, 9 formato \variable", 9 formato por defecto, 9 formato prede nido, 9 fplot2d, 45 fplot3d, 47 fprintf, 32 Frobenius, 19 fsolve, 53 funcion (programa), 22, 24 funcion exponencial, 7 funciones, 6 funciones elementales, 6 funciones matematicas, 6 function, 24 gc, 65 Getf, 25 getf, 25, 27 GIF, 49 Gomez, 2 gradiente, 55, 64

 INDICE ANAL ITICO gradiente conjugado, 65 gradiente numerico, 33 gra cas, 43 en dos dimensiones, 43 en tres dimensiones, 47 guion (programa), 22 Help, 8 help, 8 Help Dialog, 8 horner, 11 i, 9 if, 28 igual, 27 igualdades, 62 imag, 9 incremento, 29 INRIA, 1 int, 7, 8 integracion numerica, 54 interpl, 56 interpolacion, 56 linea, 56 intg, 54 inv, 16 inversa de una matriz, 16 inversa generalizada, 20 inverso aditivo, 6 isdef, 41 jacobiana, 53 jacobiana, matriz, 55 Karmarkar, 68 karmarkar, 68 leastsq, 66 leg, 45 libro, 2 l mite inferior, 29 superior, 29 linear matrix inequalities, 68 linpro, 60 LMI, 68 lmisolver, 68 lmitool, 68 log, 7 log10, 7 log2, 7 logaritmo decimal, 7 logaritmo en base dos, 7 logaritmo natural, 7 lu, 20 LU, factorizacion, 20

71

Matlab, 2 matrices, 13 matriz aleatoria, 14 de ceros, 14 de derivadas parciales., 67 de unos, 14 de nida positiva, 19 diagonal, 17 dispersa, 51 escalonada reducida por las, 17, 34 identidad, 14 inversa, 16 jacobiana, 53, 55 simetrica, 19 \vac a", 16 max, 7, 17 maximo, 7, 17 maximo comun divisor, 31 mayor, 27 mayor o igual, 27 mayusculas, 4 mean, 18 median, 18 mediana, 18 menu, 8 menor, 27 menor o igual, 27 metodos numericos, 50 min, 7, 18 m nimo, 7, 18

72 m nimos cuadrados no lineales, 66, 67 m nimos cuadrados, 54 minusculas, 4 modulo, 7, 8 modulo de un complejo, 9 multiplicacion, 6, 9, 10 de matrices, 15 elemento a elemento, 15 por escalar, 15 nd, 65 negacion, 27 no, 27 no diferenciable, 65 nombres de variables, 5 norm, 18 norma, 18 norma matricial, 18 notacion cient ca, 4, 9 numero de columnas, 19 de las, 19 numero aleatorio, 7 numero de argumentos, 35 numero de parametros, 35 o, 27 ode, 58 ones, 14 operaciones, 6 operadores logicos, 27 optim, 64 optimizacion, 60 cuadratica, 63 lineal, 60 no lineal, 64 no restringida, 64 semide nida, 68 ordenamiento, 17 parametro de entrada, 24 de salida, 24

INDICE ANAL ITICO parentesis, 6 parte imaginaria, 9 real, 9 parte entera inferior, 7, 8 superior, 7, 8 parte triangular inferior, 19 superior, 19 pegar matrices, 16 pinv, 20 plot2d, 43 plot3d, 47 polinomios, 10 poly, 10, 11 polynomial, 8 Postscrip-Latex, 49 Postscript, 49 Postscript No Preamble, 49 potencia, 6, 10, 15, 16 PPM, 49 precedencia, v ease prioridad print, 32 printf, 32 prioridad, 6 prod, 18 producto de elementos de ..., 18 producto de matrices, v ease multiplicacion ... programacion ..., v ease optimizacion ... programacion lineal, v ease optimizacion lineal programas, 22 promedio, 18 pwd, 26 qn, 65 QR, 20 qr, 20 quapro, 63

 INDICE ANAL ITICO r, 10 ra z de una matriz, 16 raiz cuadrada, 7 rand, 7, 14 rango, 16 rank, 16 real, 9 reales, 4 redondeo, 7, 8 residuo entero, 7, 8 resta, 6, 9, 10 de matrices, 15 restricciones de caja, 61, 62, 66, 67 return, 31, 32 roots, 10 Rosenbrock, 64 round, 7, 8 rref, 17 sce, 22 sci, 24 script, v ease guion select, 30 semidef, 68 seno, 7 seno hiperbolico, 7 seudoinversa, 20 seudosolucion, 23 sin, 7 sinh, 7 sistema de ecuaciones diferenciales, 58 sistema de ecuaciones no lineales, 53 sistemas de ecuaciones lineales, 50 size, 19 solucion de sistemas de ecuaciones, 20 solucion por m nimos cuadrados, 21, 23 sort, 17 sparse, 51 sparse, 51 spec, 19 spline, 56, 57 sqrt, 7 st deviation, 18 stacksize, 48 submatriz, 14 sum, 18 suma, 6, 9, 10 de elementos de ..., 18 de matrices, 15 suprimir una columna, 16 una componente, 16 una la, 16 SVD, 20 svd, 20 tama~ no de una matriz, 19 tan, 7 tangente, 7 tangente hiperbolica, 7 tanh, 7 then, 28, 30 tiempo, 50 transpuesta de una matriz, 15 tril, 19 triu, 19 valor absoluto, 6, 9 valores propios, 19 variable, 3 variables globales, 22 locales, 22 vector columna, 13, 15 vector la, 15 vectores, 13 vectores propios , 19 while, 29 who, 23 who, 4 Windows, 2 write, 32 www, 1 Xfig, 49 xtitle, 45

73

74 y, 27 zeros, 14

INDICE ANAL ITICO