Transcript
Matem´atica Num´erica Patricia Saavedra Barrera∗
1
Introducci´ on
La matem´atica num´erica es la rama de las matem´aticas que se ocupa de desarrollar m´etodos y herramientas para resolver problemas matem´aticos en forma num´erica con el fin de que los resultados cuantitativos proporcionen informaci´on sobre el comportamiento de un fen´omeno natural o de un proceso industrial, econ´omico o social. El objetivo es desarrollar algoritmos eficientes y precisos, que puedan ser implementados en un instrumento de c´alculo, cuyo desempe˜ no pueda ser analizado en detalle y cuya aplicabilidad sea lo m´as amplia posible. Los problemas que se abordan provienen de casi todas las ramas de la matem´atica tanto pura como aplicada. Hay aspectos num´ericos en la teor´ıa de n´ umeros, la combinatoria, el ´algebra abstracta y lineal; en el an´alisis, en las ecuaciones diferenciales ordinarias y parciales, en la probabilidad y en la estad´ıstica, por citar algunos ejemplos. El desarrollo vertiginoso de las computadoras trajo consigo una demanda creciente por procedimientos m´as potentes para resolver problemas cada vez m´as complejos como la predicci´on del clima, el estudio de las prote´ınas, el desarrollo de tumores cancerosos o el control de vuelo de una nave espacial. Todos estos problemas requieren de resultados cuantitativos que permitan tomar acciones, algunas a muy corto plazo, lo que implica tener m´etodos eficientes que den resultados con la precisi´on deseada en el menor tiempo posible, incluso en lo que se ha dado por llamar en tiempo real. Este campo tiene inter´es tanto para los que requieren de c´alculos como para aquellos que procesan grandes cantidades de informaci´on en forma de im´agenes, sonido o texto, ya que detr´as de cada imagen enviada y cada acceso a un buscador en la red, hay uno o varios algoritmos que requieren de miles de operaciones aritm´eticas. Desde Bacon y Galileo en el siglo XVII, el desarrollo de las ciencias naturales se ha enriquecido al combinar las metodolog´ıas experimentales con las te´oricas, estas ∗
Departmento de Matem´aticas. M´exico
Universidad Aut´onoma Metropolitana.
1
09340, Iztapalapa,
u ´ltimas muy relacionadas con la matem´atica. A partir del desarrollo de las computadoras, una tercera v´ıa que se conoce como simulaci´on num´erica ha venido a enriquecer ´ el m´etodo cient´ıfico para la generaci´on de nuevos conocimientos. Esta funciona como un laboratorio cuya herramienta principal es la computadora y en donde se puede dise˜ nar con toda precisi´on un experimento que proporciona un mejor entendimiento del fen´omeno o proceso que se estudia. Esta metodolog´ıa puede considerarse parte de la matem´atica num´erica y est´a estrechamente relacionada con varias ´areas de la computaci´on y de otras disciplinas; por ello, se habla de f´ısica, qu´ımica, biolog´ıa y hasta finanzas computacionales que se ocupan del estudio de modelos matem´aticos y algoritmos computacionales precisos y eficientes para la resoluci´on de problemas que surgen en estas ramas de la ciencia. La influencia que ha tenido la computadora en el desarrollo de la ciencia y de nuevas tecnolog´ıas corrobora la afirmaci´on visionaria que hizo John Von Neumann en 1950: ”La ciencia fue revolucionada por la invenci´on del c´alculo. El impacto de las computadoras en la ciencia ser´a al menos tan grande como el del c´alculo.”
2
Breve historia de la matem´ atica num´ erica
La matem´atica num´erica implica c´alculo con n´ umeros y registros num´ericos existen desde las primeras civilizaciones. Se han encontrado textos relacionados con problemas matem´aticos en tablillas de arcilla que provienen de excavaciones en las regiones en donde se asentaron los sumerios, caldeos, asirios y babil´onicos. Los m´as antiguos datan del u ´ltimo periodo sumerio, alrededor de 2100 A. de C., mientras que la mayor´ıa son posteriores a 600 A. de C. Las tablillas muestran que representaban los n´ umeros en un sistema posicional mixto en base 10 y 60. Con este sistema sumaban y multiplicaban con la misma facilidad con la que se hace en el sistema decimal. Manejaban las cuatro operaciones b´asicas. Para dividir se ayudaban con tablas del inverso de un n´ umero que les permit´ıa obtener el cociente aproximado de dos n´ umeros por medio de una multiplicaci´on. Hay tambi´en registro de lo que, posiblemente, es el primer algoritmo para aproximar la ra´ız cuadrada de un n´ umero. Los problemas de geometr´ıa en los que se interesaron se relacionaban con la medici´on de per´ımetros, ´areas y vol´ umenes de algunos s´olidos. Para el c´alculo del per´ımetro de un c´ırculo multiplicaban su di´ametro por 3, lo que equivale aproximar al n´ umero π por 3, ver [3]. Sobre los conocimientos matem´aticos de la cultura egipcia los primeros registros que se tienen son el Papiro de Mosc´ u y de Rhind, escritos aproximadamente en 1850 A. de C. y 1650 A. de C., respectivamente. Ambos documentos incluyen ejemplos de c´alculos que implican el manejo de ecuaciones lineales con una y dos inc´ognitas. En cuanto a geometr´ıa, determinaron con ´exito el ´area y el volumen de diversas figuras 2
geom´etricas: entre ellas, el volumen de una pir´amide truncada. Contaban con una f´ormula aproximada para calcular el ´area de un c´ırculo en el que el n´ umero π se aproxima por 3 16 . Es con los griegos que aparece por primera vez, en el siglo V A.de C., la matem´atica como una ciencia formal que utiliza el m´etodo deductivo como herramienta fundamental para probar que un resultado es verdadero. Los teoremas que se atribuyen a Tales de Mileto y a Pit´agoras son algunos de los primeros ejemplos en el que se aplica esta metodolog´ıa y el libro de los Elementos de Euclides es la muestra m´as acabada de ello. La matem´atica griega se ocup´o de estudiar problemas de geometr´ıa, aritm´etica, teor´ıa de n´ umeros y ´algebra. Arqu´ımedes (287-212) es posiblemente el m´as brillante matem´atico aplicado de la antig¨ uedad; a ´el se deben diversos ejemplos de m´etodos muy ingeniosos para aproximar la soluci´on de problemas de la hidrodin´amica y est´atica. Escribi´o m´as de diez obras; entre ellas el estudio de las c´onicas, de la esfera y el cilindro, sobre la espiral y sobre el c´alculo de vol´ umenes de revoluci´on de elipses, par´abolas e hip´erbolas que giran alrededor de un eje de simetr´ıa y el c´alculo del ´area de una par´abola. En su estudio del c´ırculo mejora sustancialmente la aproximaci´on al n´ umero π con el 1 siguiente resultado: 3 10 < π < 3 . 71 7 Desde la aparici´on de los primeros asentamientos humanos se usaron piedras, cordeles, varillas como instrumentos de medici´on y c´alculo; posteriormente, la aparici´on del ´abaco facilit´o enormemente el c´alculo aritm´etico. La introducci´on de la notaci´on decimal y del cero en el siglo XII por los ´arabes uniform´o la notaci´on y el c´alculo num´erico entre los diversos pueblos. Pero no fue sino hasta inicios del siglo XVII, con el invento de las tablas de logaritmos, que el c´alculo num´erico experiment´o un avance fundamental para el desarrollo de la ciencia en el siglo XVII y XVIII. La tabla de logaritmos fue creada por John Napier (1550-1617), quien fue de los primeros en estudiar la notaci´on y el manejo de distintas bases num´ericas como la binaria o la exponencial. Sus estudios sobre la relaci´on entre las progresiones geom´etricas y aritm´eticas lo llevaron a definir el logaritmo base 1e que public´o en 1615. Al morir Napier, su amigo Henry Briggs (1561-1630), el primer profesor saviliano de geometr´ıa de Oxford, retoma y simplifica su trabajo con el fin de generalizarlo a cualquier base, en particular la decimal, que publica en 1617 con una precisi´on de catorce cifras. Con la tablas de logaritmos, las multiplicaciones de varias cifras se convierten en simples sumas y el elevar un n´ umero a un exponente se reduce a una simple multiplicaci´on. Esto permiti´o abordar problemas de astronom´ıa, de bal´ıstica y de ´ındole pr´actico que parec´ıan incalculables. Poco despu´es, hizo su aparici´on en 1624 la regla de c´alculo que permaneci´o como instrumento personal de c´alculo por excelencia hasta la aparici´on de las primeras calculadoras de bolsillo en 1970. Con estos antecedentes los precursores del c´alculo: Descartes, Fermat, Pascal, Wallis, Barrow y Gregory, entre otros, desarrollaron la geometr´ıa anal´ıtica, diversos 3
Figura 1: John Napier m´etodos para determinar la ecuaci´on de la tangente a una curva dada; el c´alculo de m´aximos y m´ınimos y el c´alculo de cuadraturas, entre otros temas. Estos resultados prepararon el terreno para que, a fines del siglo XVII, Leibnitz y Newton los sistematizaran, relacionaran y unificaran en una metodolog´ıa poderos´ısima: el c´alculo diferencial e integral. Isaac Newton (1642-1727) naci´o en Woolsthorpe, Inglaterra. Sus aportaciones son fundamentales para definir con precisi´on el concepto de velocidad inst´antanea, la fuerza de atracci´on entre los cuerpos, la fuerza gravitacional y la descomposici´on de la luz, entre otras muchas cosas. Desde el punto de vista num´erico, Newton invent´o varios algoritmos que se siguen utilizando hasta nuestros d´ıas, como el m´etodo para encontrar los ceros de una funci´on o la aproximaci´on de derivadas por medio de diferencias divididas o la construcci´on de polinomios de interpolaci´on. En esa misma ´epoca se introdujo el uso de series de potencias que permiti´o aproximar localmente con la precisi´on deseada a un gran conjunto de funciones y abri´o la posibilidad de atacar problemas de mayor complejidad que aparecen en el estudio de fen´omenos f´ısicos como las vibraciones de una cuerda, la difusi´on de calor o el comportamiento de los fluidos. Por ejemplo, James Gregory (1638-1675) en 1671 determin´o la serie de potencias para la funci´on Arctang(x) lo que permiti´o aproximar,
4
Figura 2: Isaac Newton tanto como se deseara, el valor de π π 1 1 1 = 1 − + − + .... 4 3 5 7 Siguiendo el ejemplo de Newton, varios matem´aticos de primer orden contribuyeron a la soluci´on de problemas num´ericos: Leonard Euler (1707-1783) utiliz´o la serie de Taylor para crear el m´etodo que lleva su nombre y que aproxima la soluci´on de problemas con condiciones iniciales en ecuaciones diferenciales ordinarias; Joseph Louis Lagrange (1736-1813) di´o una forma de construir polinomios de interpolaci´on. Karl Friedrich Gauss (1777-1855) contribuy´o de manera notable con el m´etodo de m´ınimos cuadrados y sus f´ormulas de cuadratura que abren la puerta a la utilizaci´on de familias de polinomios ortogonales para la aproximaci´on de integrales; sus aportaciones al estudio de la soluci´on de sistemas de ecuaciones lineales siguen siendo fundamentales en la ense˜ nanza del ´algebra lineal. Tambi´en se le atribuye la invenci´on de la transformada r´apida de Fourier, algoritmo que pas´o desapercibido hasta su redescubrimiento por Cooley y Tukey en 1963. Durante la segunda mitad del siglo XIX, la matem´atica entr´o en una etapa de rigorismo necesario para cimentar con firmeza su amplio cuerpo de conocimientos. 5
Weierstrass y Cauchy entre otros, son los responsables del concepto riguroso de l´ımite que es fundamental en el an´alisis. Esto explica que la obtenci´on de c´alculos num´ericos pasara a segundo t´ermino aunque el desarrollo del an´alisis funcional est´a estrechamente ligado a la aproximaci´on de funciones en espacios normados. Por u ´ltimo, hay que mencionar que, a principios del siglo XX, Richardson aplic´o por primer vez el m´etodo de diferencias finitas para el an´alisis de esfuerzos en distintos materiales y que, en 1928 apareci´o el famoso art´ıculo de Courant, Friedrichs y Lewy sobre el an´alisis de estabilidad de estos m´etodos. Para saber m´as sobre la historia del an´alisis num´erico, consultar [8].
Desarrollo de las computadoras En contraste con el desarrollo acelerado de la ciencia en el siglo XIX y la primera mitad del siglo XX, los instrumentos de c´alculo permanecieron sin modificaci´on alguna hasta mediados del siglo XX. Charles Babbagge en el siglo XIX fue el primero en tener la idea de dise˜ nar un m´aquina que realizara c´alculos autom´aticos, pero sus ideas no fueron comprendidas por sus contempor´aneos y tuvieron que esperar hasta inicios de la Segunda Guerra Mundial para dar frutos. La primera computadora, llamada Z3, apareci´o en Alemania en 1939 y fue constru´ıda por el ingeniero Konrad Zuse; usaba dispositivos electromec´anicos y representaba a los n´ umeros en binario en notaci´on de punto flotante, t´ermino que ser´a definido m´as adelante. Tambi´en durante esos a˜ nos, los brit´anicos construyeron una gran computadora para descifrar los c´odigos de encriptamiento de la m´aquina alemana Enigma.
Figura 3: La computadora ENIAC La primera computadora de uso general fue la ENIAC, construida por Prosper Eckert y John Mauchly en la Universidad de Pennsylvania de 1943 a 1945. Aunque esta computadora era electr´onica, cada vez que se usaba en un problema nuevo hab´ıa 6
que reprogramarla lo que implicaba reconectar distintos circuitos el´ectricos. En el verano de 1944, Von Neumann visit´o la construcci´on de la ENIAC como asesor de la marina. El problema de proporcionar a la computadora las instrucciones a seguir sin necesidad de desconectar los circuitos despert´o su inter´es. En 1945 Von Neumann escribi´o un reporte sobre el dise˜ no l´ogico de una computadora nueva, llamada EDVAC; en ´el plasma las ideas que hab´ıan surgido en sus discusiones con Eckhert, Mauchly y Burks, entre otros, ver [9], y dibuja un diagrama de lo que debe ser la estructura de una computadora: una unidad de c´alculo, una unidad central que organice las tareas que lleva a cabo la computadora y una memoria en la que se almacenen los programas. Este documento tendr´a una influencia determinante en el desarrollo posterior de las computadoras; su concepci´on permaneci´o vigente hasta la aparici´on de las computadoras en paralelo. Al mismo tiempo Von Neumann se interes´o en la confiabilidad de los resultados que se obtienen en una computadora, cuesti´on nada evidente en esos a˜ nos. Estudi´o la propagaci´on de errores en los c´alculos num´ericos y el dise˜ no de algoritmos menos sensibles a ´estos que son aspectos fundamentales de la matem´atica num´erica.
Figura 4: John Von Neumann Antes de la aparici´on de las computadoras, a la matem´atica num´erica se le consideraba parte del an´alisis ya que se ocupaba principalmente del estudio de la aproximaci´on num´erica de la soluci´on de problemas continuos. A partir de 1945, el desarrollo de la matem´atica num´erica se aceler´o, a la par que lo hizo la computadora, 7
y desbord´o su marco inicial para ocuparse de todo tipo de problemas cuyo estudio y soluci´on pod´ıan beneficiarse del uso intensivo de esta herramienta. En cuanto a las publicaciones, Teor´ıa Matem´atica de la Computaci´on de la Sociedad Matem´atica Americana inici´o su publicaci´on en 1943, pero no fue sino hasta la d´ecada de los cincuenta que aparecieron las primeras revistas especializadas totalmente dedicadas a este nuevo campo; entre las primeras est´an el SIAM Journal en 1953 y el Numerische Mathematik en 1959. De esta forma el n´ umero de matem´aticos, ingenieros y cient´ıficos que se interesaron en crear algoritmos que aprovecharan la capacidad de las computadoras se increment´o tanto que, a principios de los sesenta, la matem´atica num´erica se hab´ıa convertido en una rama independiente de la matem´atica. Actualmente, hay m´as de una docena de revistas especializadas en este tema, sin contar las que se publican en otras disciplinas. La velocidad de las computadoras se ha incrementado en un factor de 109 desde su creaci´on, lo que se ha convertido en un reto constante para generar nuevos algoritmos que puedan aprovechar y explotar dicho potencial. Hoy en d´ıa, se puede afirmar que, gracias al esfuerzo de muchos investigadores, buena parte de los problemas cl´asicos de la f´ısica pueden ser resueltos computacionalmente en forma confiable y eficiente; pero a´ un restan otros problemas, como el estudio de la turbulencia o del clima, que junto con otros nuevos, como el dise˜ no de prote´ınas sint´eticas o de nuevos materiales siguen siendo un reto para el ingenio de las nuevas generaciones de analistas num´ericos.
3
La matem´ atica num´ erica en M´ exico
La matem´atica num´erica en M´exico tiene las mismos antecedentes que cualquier otro campo de la matem´atica en M´exico. Una de las caracter´ısticas del desarrollo de esta disciplina en nuestro pa´ıs es su falta de continuidad, por lo que se puede hablar de tres periodos ajenos: el precolombino, en el que destaca la matem´atica maya, el colonial, consultar el libro de E. Trabulse, ver [19], y el posrevolucionario a partir de 1930 a nuestros d´ıas, en el que se profesionaliza la actividad matem´atica en M´exico. A continuaci´on se hace una breve presentaci´on de esta u ´ltima etapa. Despu´es de la Revoluci´on, Don Sotero Prieto (1884-1935), profesor de la Escuela Nacional Preparatoria y de la Escuela Nacional de Ingenieros, reuni´o a un grupo de j´ovenes interesados en la matem´atica para profundizar ciertos temas e intercambiar ideas. La importancia de este grupo es que varios de sus miembros jugaron un papel determinante en el impulso de la actividad matem´atica en M´exico. Uno de estos j´ovenes era Alfonso N´apoles G´andara, quien en 1930 recibi´o una beca para realizar estudios de matem´aticas en el M.I.T. En 1935 surgieron las carreras de matem´atico y f´ısica, dentro del Departamento de Ciencias F´ısicas y Matem´aticas de la UNAM y, en noviembre de 1938, se aprob´o la creaci´on de la Facultad de Ciencias, el Instituto 8
de F´ısica y el Instituto de Matem´aticas. Estas instituciones junto con la Sociedad Matem´atica Mexicana creada en 1943, lograron contagiar a profesores y alumnos de distintos lugares de nuestro pa´ıs del gusto por las distintas ramas de las matem´aticas y ayudaron a crear, antes de 1970, la carreras de matem´aticas en varias universidades de provincia como Nuevo Le´on, Puebla, Veracruz y Yucat´an.
Figura 5: N´apoles G´andara N´apoles G´andara fue el director de Instituto de Matem´aticas desde su inicio hasta 1965. Organiz´o el Instituto en tres ramas generales: Matem´atica Pura, dirigida por Alberto Barajas y Roberto V´azquez; L´ogica y Fundamentos, dirigida por Francisco Zubieta, y Matem´atica Aplicada cuyo responsable era Carlos Graef. En los siguientes veinte a˜ nos, s´olo la primera rama lograr´ıa consolidarse con el apoyo decidido de Solomon Lefschetz, investigador de la Universidad de Princeton. Sus visitas a M´exico, a partir de 1945 hasta 1966, influyeron considerablemente en la formaci´on de investigadores mexicanos en Estados Unidos en tres ´areas relevantes de la matem´atica: geometr´ıa algebraica, ecuaciones diferenciales y topolog´ıa algebraica. El desarrollo de la matem´atica aplicada en M´exico fue m´as lento y m´as disperso ya que esta actividad no es exclusiva de los matem´aticos, sino que tambi´en se re9
troalimenta con el trabajo realizado por los f´ısicos, geof´ısicos e ingenieros. En consecuencia, otras instituciones influyeron en el avance de la matem´atica aplicada: la Facultad de Ciencias (FC), el Instituto de Investigaciones en Matem´aticas Aplicadas y Sistemas (IIMAS), el Instituto de Ingenier´ıa y el Instituto de Geof´ısica (IG) de la UNAM. En este instituto hay que destacar el trabajo del Dr. Ismael Herrera, matem´atico egresado de la Facultad de Ciencias y uno de los primeros en obtener un doctorado en matem´aticas aplicadas. El Dr. Herrera contribuy´o en el desarrollo de los m´etodos variacionales en la soluci´on de ecuaciones en derivadas parciales y su aplicaci´on a la ingenier´ıa. Tambi´en hay que mencionar la importancia de la labor realizada por el Dr. Julian Adem y el Dr. Arturo Rosenblueth de la UNAM y el trabajo pionero realizado por el Ing. Manuel Cerrilla quien public´o desde 1945 varios trabajos sobre algunos m´etodos matem´aticos aplicados a problemas el´ectricos, de electromagnetismo y procesamiento de se˜ nales. Vale la pena recordar el desarrollo del IIMAS ya que su historia est´a entrelazada con el desarrollo de la computaci´on en M´exico. El Centro de C´alculo Electr´onico (CCE), fue fundado en 1958, a˜ no en que se instala la primera computadora en la UNAM y en el pa´ıs, con el fin de utilizarla para el avance de la ciencia en M´exico. A partir de la creaci´on del CCE, el n´ umero de cient´ıficos y profesionales que usaron esta nueva herramienta como parte de sus investigaciones no dej´o de crecer y mostr´o la necesidad de formar recursos humanos en computaci´on y ´areas afines en el extranjero. En 1967 el Centro se moderniza, adquiere una computadora con tecnolog´ıa muy avanzada para su tiempo y su uso se difunde r´apidamente, pasando de 60 a 2000 usuarios activos. El CCE cambia de nombre (CIMASS) y funciones en 1970, reconociendo entre sus actividades la investigaci´on en c´omputo y en estad´ıstica. En 1973 se divide el CIMASS para separar la actividad de investigaci´on de la de servicios y se crea el CIMAS que, posteriormente en 1976, se convertir´a en el IIMAS. Este instituto ha desempe˜ nado un papel fundamental en el desarrollo de la matem´atica aplicada y la estad´ıstica en M´exico. En 1962 el Centro de Estudios Avanzados del Instituto Polit´ecnico Nacional (IPN) inici´o sus actividades. El Dr. Jos´e Adem, hasta entonces investigador del Instituto de Matem´aticas de la UNAM, se encarg´o de dirigir este nuevo centro que ha tenido gran impacto en el desarrollo matem´atico de nuestro pa´ıs. Debe mencionarse que un poco antes, en 1961, se estableci´o la Escuela Superior de F´ısica y Matem´aticas en el IPN, cuyos primeros profesores eran egresados de la Facultad de Ciencias de la UNAM. En 1974 se crea la Universidad Aut´onoma Metropolitana, cuyo Departamento de Matem´aticas en la Unidad Iztapalapa, ha formado a un buen n´ umero de matem´aticos aplicados en el pa´ıs en las ´areas de control, f´ısico matem´aticas, sistemas din´amicos y an´alisis num´erico. En 1980 se funda el Cimat en Guanajuato, instituci´on que ha sido un foco de desarrollo y de atracci´on de los matem´aticos en provincia. En 1991 se llev´o a cabo, a partir de una iniciativa entusiasta de varios investi10
gadores de la Facultad de Ciencias, entre ellos los doctores Pablo Barrera y Jes´ us L´opez Estrada, y el Dr. Humberto Madrid de la Universidad Aut´onoma de Coahuila, la primera Escuela de Optimizaci´on y An´alisis Num´erico con el fin de difundir y actualizar entre los profesores e investigadores de las universidades de provincia, los avances en estas ´areas. Hasta ahora, se han llevado a cabo diecisiete escuelas lo que ha contribuido a la consolidaci´on de la matem´atica aplicada en nuestro pa´ıs. Tambi´en hay que se˜ nalar en esta direcci´on el trabajo de divulgaci´on y promoci´on realizado por el Dr. Diego Bricio Hern´andez, ver [11]. Actualmente, hay alrededor de doscientos investigadores activos en el ´area de matem´aticas aplicadas y, entre ellos, cerca de cincuenta son en matem´atica num´erica, sin contar aquellos que provienen de otras disciplinas. Las ´areas que se cultivan son la optimizaci´on, la soluci´on num´erica de ecuaciones diferenciales ordinarias y parciales, el ´algebra lineal num´erica y los m´etodos num´ericos aplicados a la probabilidad y a la estad´ıstica. Se encuentran grupos consolidados en el CIMAT; en la Facultad de Ciencias, en el IIMAS, en los Institutos de Ciencias de la Atm´osfera, Geof´ısica e Ingenier´ıa de la UNAM; en el ITAM; en el Departamento de Matem´aticas de la Unidad Iztapalapa de la UAM y en las Facultades de Ciencias F´ısico-Matem´aticas de la BUAP, de la U. M. de S. N. Hidalgo y de la U.A.N.L.
4
Matem´ atica num´ erica y las computadoras
Hay muchas formas de representar los n´ umeros, desde la antig¨ uedad se conocen las ventajas de utilizar potencias de 10 lo que da un sistema posicional decimal en el que es importante la posici´on que ocupa un d´ıgito. Usar las potencias de diez es natural porque los humanos tienen diez dedos en la mano, pero hay otros sistemas que utilizan las potencias de otros n´ umeros. Los babil´onicos ten´ıan un sistema sexagesimal y quedan vestigios de ello en nuestra forma de medir el tiempo. Para las computadoras lo m´as natural es usar potencias de dos, ya que almacenan los caracteres en bytes que a su vez se componen de 8 bits. Cada bit puede tomar dos valores, uno si est´a prendido y cero si est´a apagado. Un byte permite representar 256 caracteres distintos. Una palabra de computadora se compone de 4 bytes y una palabra doble de 8 bytes consecutivos. Al usar potencias de dos se dice que se trabaja en base binaria en lugar de la base decimal. Algunas computadoras modernas trabajan en base octal, base 8, o hexadecimal que es en base 16. La manera m´as f´acil de representar un n´ umero real en una computadora es en notaci´on de punto flotante. A continuaci´on, por simplicidad, se presenta esta notaci´on en base decimal aunque las computadoras usan otras bases. La notaci´on es la siguiente: n x=+ − .d1 d2 d3 ..........dk dk+1 ...... × 10 , 11
donde n es un n´ umero entero, dk es un n´ umero entero entre 0 y 9, y d1 es distinto de cero. La sucesi´on de d´ıgitos dk se le conoce como la mantisa y al entero n como el exponente. Como la representaci´on de un n´ umero real puede tener mantisa infinita y n puede tomar cualquier valor en los enteros, es imposible que un instrumento de c´alculo pueda representar cualquier n´ umero en forma exacta; cada instrumento de c´alculo tiene una cota superior e inferior para n y s´olo puede almacenar un n´ umero finito de d´ıgitos k de la mantisa. Si se almacena un real con exponente mayor a n, la computadora marca un error que se conoce como overflow o desbordamiento superior. Si el real tiene como exponente un entero menor a la cota inferior, la m´aquina genera un error de underflow o desbordamiento inferior. Cuando la mantisa tiene un n´ umero de d´ıgitos mayor a k entonces la computadora redondea el n´ umero y almacena s´olo k d´ıgitos. El procedimiento de redondear consiste en lo siguiente: primero se observa el valor del d´ıgito dk+1 ; si ´este es menor o igual a 4 se almacenan u ´nicamente los primeros k d´ıgitos; si dk+1 es mayor o igual a 5 se suma un uno al k ´esimo d´ıgito y se almacenan los k d´ıgitos del n´ umero resultante. Para representar un n´ umero en punto flotante se puede usar una palabra simple o doble, la diferencia estriba en el n´ umero de d´ıgitos k de la mantisa que almacenan. Si se usan palabras simples para representar a los reales se dice que se calcula en simple precisi´on o en doble precisi´on en caso de utilizar palabras dobles. ¿De qu´e tama˜ no es el error que se comete al representar cualquier n´ umero real en una computadora? Para estimarlo den´otese por f l(x) la representaci´on en punto flotante de un n´ umero real x en una computadora con una mantisa de k d´ıgitos. Se usar´a el error relativo Er para estimar el error y se puede probar que Er =
|x − f l(x)| ≤ .5 × 101−k . |x|
A este error se le llama error de redondeo. En consecuencia, cuando se resuelve un problema con una computadora, la mayor parte de las veces, se resuelve un problema aproximado y no el original; los resultados obtenidos tendr´an, en el mejor de los casos, la precisi´on de la computadora.
Aritm´ etica en punto flotante ¿Qu´e sucede cuando realizamos las cuatro operaciones de la aritm´etica con f l(x) en lugar de x? Uno de los problemas que se presentan es que las cuatro operaciones no son cerradas en el conjunto de f l(x) con x real. Por ejemplo, en caso de que haya desbordamiento, como cuando se multiplican dos n´ umeros muy grandes, no podr´a llevarse a cabo la operaci´on y la computadora marcar´a el error de desbordamiento. Aunque pueda llevarla a cabo, el resultado de la operaci´on puede tener mantisa de m´as 12
de k d´ıgitos por lo que hay que redondearla. El resultado de las cuatro operaciones en notaci´on punto flotante es de la forma: x+ ≈ f l(f l(x)+ − f l(y)), −y x × y ≈ f l(f l(x) × f l(y)), x/y ≈ f l(f l(x)/f l(y)). A primera vista, parecer´ıa que los errores de redondeo se ir´an acumulando a medida que se realicen las operaciones. Si as´ı fuera, ser´ıa imposible usar una computadora o cualquier otro instrumento de c´alculo para aproximar la soluci´on de un problema que involucre un gran n´ umero de operaciones. Por fortuna los errores de redondeo se comportan en forma estad´ıstica y no determinista; tienen un distribuci´on normal por lo que, en la mayor´ıa de los c´alculos, el resultado conserva la misma precisi´on que la de la computadora. Sin embargo, en ocasiones, el error de redondeo crece sin parar, deteriorando la precisi´on de los resultados. Numerosos estudios realizados por Von Neumann, Wilkinson, Forsythe y Henrici entre otros, documentaron ampliamente este fen´onemo al que llamaron inestabilidad num´erica, ver [6]. Los estudios muestran que este comportamiento no s´olo es resultado del n´ umero de d´ıgitos de precisi´on de la computadora, sino tambi´en de las caracter´ısticas del problema mismo y del algoritmo que se utilice. Uno de los aspectos que debe tomar en cuenta, cualquier persona interesada en aproximar num´ericamente la soluci´on de un problema, es qu´e tan sensible es ´este y el algoritmo que va a utilizar al error de redondeo. La notaci´on en punto flotante era ya est´andar a mediados de la d´ecada de 1950, pero cada fabricante desarrollaba su propio sistema y, en ocasiones, distintos programas en una misma computadora lo hac´ıan en forma diferente, con resultados desastrosos. En 1985 se aprob´o un acuerdo que se conoce por Estandar IEEE, tanto para la notaci´on, como para la aritm´etica en punto flotante en binario en todas las computadoras que se produc´ıan en Estados Unidos. En 1987 se adopt´o ese mismo est´andar para cualquier base, ya que las calculadoras utilizan base decimal, y en 1989 se adopt´o a nivel internacional. A pesar de la aprobaci´on del est´andar IEEE la falta de consistencia entre los distintos softwares y el uso de programas anteriores a la aprobaci´on del acuerdo, han ilustrado recientemente los problemas que provoca un manejo descuidado de la aritm´etica en punto flotante. Entre los m´as graves est´an, en 1996 la autodestrucci´on del cohete Arianne V; fallas en el sistema de antimisiles Patriot durante la Guerra del Golfo y en 1997, un mal manejo del error de redondeo a nivel de hardware, oblig´o a Intel a sacar del mercado sus microprocesadores Pentium Pro y Pentium II, con un desplome inmediato en la cotizaci´on de sus acciones, ver [16] .
13
Figura 6: El cohete Arianne V
Problemas bien planteados No todos los problemas matem´aticos pueden ser resueltos con la ayuda de una computadora. El matem´atico Jacques Hadamar, a principios del siglo XX, precis´o las caracter´ısticas que debe tener un problema para estar bien planteado: la soluci´on debe existir, aunque sea localmente, ser u ´nica, y la soluci´on debe depender en forma continua de los datos. Si un problema no tiene soluci´on, no tiene sentido aproximarla num´ericamente. Por ejemplo, en el caso que se desee calcular una integral divergente. Si la soluci´on no es u ´nica, se corre el riesgo de obtener una soluci´on num´erica que no aproxime a soluci´on alguna. Por ejemplo, en el caso de que nos interese aproximar la soluci´on de una ecuaci´on diferencial de la forma y 0 (x) = f (x, y(x)) x ∈ (0, 1), y(0) = y0 . 14
Si no se puede garantizar unicidad en (0, 1), puede suceder que para algunos valores de xi la soluci´on num´erica se aproxime a una de las curvas integrales, mientras que para otros valores aproxime a otra curva integral. Al graficar los valores obtenidos, no corresponder´an a soluci´on alguna. Por u ´ltimo, la dependencia continua de la soluci´on respecto a los datos significa que si se modifican ´estos un poco, la soluci´on del problema perturbado es muy cercana a la del problema original. ¿De qu´e manera se decide matem´aticamente que dos soluciones se parecen? Este es un problema delicado que hay que precisarlo en cada caso. Para sistemas lineales de ecuaciones se tendr´a una definici´on distinta que para el caso de una ecuaci´on diferencial. Probar que un problema est´a bien planteado no es tarea f´acil y en algunos casos es todav´ıa un problema abierto a la investigaci´on. Si se sabe de antemano que un problema est´a bien planteado, es muy probable que la aproximaci´on computacional arroje resultados confiables. Sin embargo, aunque no se cumpla la tercera condici´on se puede resolver el problema computacionalmente siempre y cuando se usen los algoritmos adecuados. De hecho, mucha de la investigaci´on que se realiza hoy d´ıa est´a relacionada con el estudio de problemas mal planteados como los que aparecen con frecuencia en la soluci´on de problemas inversos.
Algoritmos A diferencia de otras ramas de la matem´atica, probar que un problema est´a bien planteado no es suficiente para el matem´atico num´erico, ya que debe encontrar una forma de aproximar, tanto como quiera, la soluci´on. Generar buenos algoritmos num´ericos es un arte. A veces estos son resultado directo de la teor´ıa, pero otras veces se requiere formular de manera distinta el problema; por ejemplo, en lugar de buscar los ceros de una funci´on se buscan los puntos fijos de un problema equivalente. As´ı como hay problemas que no son resolubles en una computadora tambi´en hay algoritmos cuya implementaci´on computacional no se recomienda debido a su poca eficiencia. La eficiencia de un algoritmo se mide por el n´ umero de operaciones num´ericas b´asicas que involucra. Los algoritmos se clasifican seg´ un el n´ umero de operaciones que requieren realizar para obtener el resultado deseado. Si un problema con n inc´ognitas requiere n operaciones se dice que es un algoritmo de complejidad lineal, si es de ns , para un natural s, se dice que es un algoritmo polinomial. Pero algunos algoritmos son de tipo exponencial o factorial, como el c´alculo de determinantes, lo que los hace incosteables; entre dos algoritmos con la misma precisi´on se prefiere aquel que requiera del menor n´ umero de operaciones. Hay algoritmos que en un paso obtienen en principio, si no hubiera error de redondeo, la soluci´on exacta. Mientras que hay otros m´etodos llamados iterativos, que a partir de una aproximaci´on inicial x0 , generan una sucesi´on de valores {xn }∞ n=1 que 15
convergen a la soluci´on cuando n tiende a infinito. Estos m´etodos son en muchos casos la u ´nica alternativa para resolver problemas no lineales. La desventaja principal de este tipo de m´etodos es que la sucesi´on debe interrumpirse antes de la convergencia, ya que cualquier algoritmo que se implemente en una computadora debe constar de un n´ umero finito de pasos. En consecuencia, s´olo se genera un n´ umero finito de elementos de la sucesi´on. Esto introduce un error en los c´alculos que se busca sea menor o igual a la precisi´on de la computadora.
Campos de la matem´ atica num´ erica El n´ umero de l´ıneas de investigaci´on que se trabajan actualmente en la matem´atica num´erica es muy grande y es imposible presentarlas todas en un texto introductorio y de car´acter general como ´este. Una posibilidad era mencionar a una buena parte de ellas de manera muy somera y con una gran lista de referencias. Para aquellos interesados en conocer una lista completa, consultar el art´ıculo de divulgaci´on de Trefethen, ver [20]. Otra opci´on era exponer unos cuantos temas de manera m´as amplia, buscando transmitir en la exposici´on los aspectos que le interesan al matem´atico num´erico. Se opt´o por esta u ´ltima. La selecci´on de temas se hizo tomando en cuenta que estos son algunos de los que m´as se han trabajado en M´exico: ´algebra lineal num´erica, soluci´on de ecuaciones diferenciales ordinarias y parciales y optimizaci´on.
5
Algebra lineal num´ erica
El ´algebra lineal num´erica es uno de los pilares que sostienen a la matem´atica num´erica y aplicada. Algunos de los temas principales de investigaci´on son la soluci´on de sistemas lineales, la aproximaci´on por m´ınimos cuadrados lineales y la obtenci´on de vectores y valores propios. S´olo se presenta la soluci´on num´erica de sistemas de ecuaciones lineales, los otros temas se pueden consultar en [7]. Una gran cantidad de aplicaciones de la matem´atica involucran la soluci´on de sistemas lineales por lo que se ha dedicado mucho trabajo y esfuerzo a la obtenci´on de m´etodos eficientes y confiables para su resoluci´on. Un sistema de ecuaciones lineales de n ecuaciones con n inc´ognitas es de la forma a11 x1 + a12 x2 + · · · + a1n xn = b1 , a21 x1 + a22 x2 + · · · + a2n xn = b2 , ... an1 x1 + an2 x2 + · · · + ann xn = bn .
16
Los sistemas lineales aparecen en el c´alculo de m´aximos y m´ınimos de una funci´on cuadr´atica de varias variables, en la soluci´on de problemas de optimizaci´on lineal con restricciones lineales, en la aproximaci´on de funciones por medio de polinomios de interpolaci´on o por m´ınimos cuadrados y en la soluci´on num´erica de problemas integrales o en derivadas parciales, entre otros. En este u ´ltimo caso implican la soluci´on de miles de sistemas de ecuaciones lineales con miles de variables.
Figura 7: Carl Friedrich Gauss Hay m´etodos directos e iterativos. Entre los primeros, el m´as general es el m´etodo de eliminaci´on de Gauss creado por el gran matem´atico Gauss. El sistema de ecuaciones anterior se puede escribir matricialmente de la forma − → → A− x = b,
(1)
− → → con A una matriz cuadrada de n renglones con n columnas y − x y b son vectores en Rn . La idea del m´etodo de Gauss es la siguiente: transformar el sistema (1) en un sistema triangular superior, ver glosario, por medio de operaciones elementales que son aquellas que no alteran la soluci´on del sistema. La soluci´on del sistema resultante − → → A0 − x = b0 es m´as f´acil de obtener que la del sistema original, ya que al resolver el sistema de abajo hacia arriba, el problema se reduce a resolver n ecuaciones con una 17
sola inc´ognita. Adem´as, el m´etodo permite comprobar que la soluci´on es u ´nica ya 0 que basta que las Aii 6= 0, para i = 1, . . . , n para que esto se cumpla. El resto de los m´etodos directos se basan en la factorizaci´on de la matriz A en matrices con estructura m´as sencilla. Por ejemplo, cuando el sistema se puede transformar a la forma triangular, sin intercambio de renglones, A puede escribirse como producto de dos matrices, es decir A = LU, con L una matriz triangular inferior, con unos en la diagonal y U es una matriz triangular superior. Esta factorizaci´on da lugar a otros m´etodos que resuelven de manera m´as eficiente los sistemas lineales que tienen una estructura especial. Por ejemplo, el m´etodo de Cholesky que se usa cuando la matriz asociada es sim´etrica y positiva definida, ver glosario. Otro aplicaci´on interesante de esta factorizaci´on aparece en la soluci´on num´erica de ecuaciones en derivadas parciales y que da lugar a la soluci´on de sistemas ralos, llamados as´ı porque un buen n´ umero de los coeficientes Aij son iguales a cero. Una de las ventajas de los m´etodos de factorizaci´on es la reducci´on del n´ umero de operaciones necesarias para obtener la soluci´on. El m´etodo de eliminaci´on de Gauss 3 requiere de aproximadamente n3 operaciones para resolver un sistema de n ecuaciones con n inc´ognitas, mientras que la factorizaci´on de Cholesky necesita la mitad. ¿Qu´e tan precisa es la soluci´on obtenida por computadora al usar m´etodos exactos? Para responder esta pregunta se requiere definir con cuidado lo que significa que dos elementos de Rn est´en cerca. La forma de hacerlo es a trav´es de una funci´on que se → le llama norma. Una norma es una funci´on ||.|| que a cada vector − x le asocia un real mayor o igual a cero que satisface ciertas P propiedades, ver [7]. La norma usual es la − euclideana que se define por ||→ x ||2 = [ ni=1 x2i ]1/2 . − → → Si − x es la soluci´on del problema original y x∗ es la soluci´on que se obtiene en la computadora, se puede probar que el error relativo es igual a − → → ||− x − x∗ ||2 ≤ C × 10−k , (2) − → || x ||2 con k la precisi´on de la computadora y C > 0 una constante que depende del producto de la norma de A y de su inversa. Este producto juega un papel muy importante y se conoce con el nombre de n´ umero de condici´on de la matriz A y se denota por cond(A). Si cond(A) es grande, tambi´en C lo es, y se dice que la matriz A est´a mal condicionada. En este caso, los resultados que se obtengan con el m´etodo de eliminaci´on de Gauss, no van a tener la precisi´on esperada ya que, por cada orden de magnitud que tenga C, se pierde un orden de precisi´on en los resultados. Esto ha motivado la creaci´on de m´etodos especiales que no sean tan sensibles a los errores de redondeo como m´etodos de precondicionamiento, para enterarse m´as sobre este tema, consultar [7] y [10]. Tambi´en se usan m´etodos iterativos para resolver los sistemas lineales de ecuaciones. Hay muchas formas de generar la sucesi´on. Una clase de m´etodos iterativos 18
se basa en la idea de transformar el problema en un problema de punto fijo de una funci´on G, el punto que bajo G va a dar a s´ı mismo. La sucesi´on se genera de la → → − → siguiente forma: dado − x0 , − x− etodo de Jacobi, Gauss-Seidel y relak+1 = G(xk ). El m´ jaci´on son algunos ejemplos. Otros m´etodos iterativos m´as recientes como gradiente conjugado o GMRES se presentan como buenas alternativas para reducir el error de redondeo y el tiempo de c´omputo, cuando el n´ umero de ecuaciones es muy grande, ver [10] para m´as informaci´on. En M´exico han hecho contribuciones a este campo: H. Madrid de la Universidad Aut´onoma de Coahuila, J. L. Morales del ITAM, M.L. Sandoval de la UAMIztapalapa, P. Barrera y J. L´opez E. de la FC y S. G´omez del IIMAS en la UNAM, entre muchos otros.
6
Ecuaciones diferenciales ordinarias
La formulaci´on matem´atica de muchos problemas de aplicaci´on se hace a trav´es de sistemas de ecuaciones diferenciales ordinarias (EDO) o parciales (EDP). Ejemplo de los primeros son los problemas de mec´anica, de circuitos el´ectricos y biol´ogicos en los que la soluci´on s´olo depende de una variable que puede ser espacial o temporal. La forma general de un problema de EDO con condiciones iniciales es y 0 (t) = f (t, y(t)) a < t < b,
y(a) = y ∗ ,
(3)
cuya soluci´on existe y es u ´nica si f es continua y Lipschitz en y en el intervalo [a, b]. Una forma de aproximar la soluci´on en un n´ umero finito de puntos del intervalo [a, b] la di´o Euler; la interpretaci´on geom´etrica es la siguiente: dado que se conoce el valor de la soluci´on y de su derivada en t0 = a, trazar la recta que pasa por (t0 , y ∗ ) con pendiente y 0 (t0 ) = f (t0 , y ∗ ), la ordenada del punto de intersecci´on de esta recta con la recta t = t1 nos da una aproximaci´on al valor de y(t1 ) y as´ı sucesivamente. Hay muchas formas de seleccionar los puntos ti , la m´as sencilla es usar puntos cuya distancia entre dos consecutivos sea la misma: h = |ti − ti−1 |. Entonces cada yi se obtiene por el esquema de Euler de la siguiente forma: dado N ∈ N , sea h = (b−a)/N y def´ınanse los siguientes puntos t0 = a, ti = a + ih, para i = 1, . . . , N y si y0 = y ∗ yi = yi−1 + hf (ti−1 , yi−1 ),
i = 1, . . . , N.
(4)
Otra manera de obtener el m´etodo de Euler es por medio de la serie de Taylor lo que permite estimar el error que se comete al aproximar y(ti ) por Euler cuando los valores anteriores y(tj ) con j = 1, . . . , i − 1 son exactos. Este error se conoce con el nombre de error local del m´etodo y es aproximadamente h para Euler. Obs´ervese que al aplicar el m´etodo, se usa el valor aproximado yi−1 para calcular el valor yi , por lo 19
que el error se va acumulando a medida que se generen m´as puntos; a este error se le conoce con el nombre de error global. Se dice que una soluci´on aproximada {yi }ni=1 converge a la soluci´on y(t) si el error global en el punto tn converge a cero cuando h tiende a cero. Es f´acil probar la convergencia del m´etodo de Euler cuando la funci´on f es suficientemente regular en el intervalo [a, b]. M´etodos m´as precisos se pueden obtener tomando m´as t´erminos de la serie de Taylor, pero tienen el inconveniente de que deben existir las derivadas de orden superior de f y, aunque existan, su c´alculo puede llegar a ser muy engorroso. Una forma de darle vuelta a esta dificultad, es aproximar las derivadas de orden superior de f por aproximaciones de la serie de Taylor en dos variables, lo que ha dado lugar a los conocidos m´etodos Runge-Kutta (RK) de distinto orden, el m´as popular es el de orden cuatro (RK4). En la Figura 8 se comparan la soluci´on exacta y las aproximadas que se obtienen con Euler y RK4 con h = 0.05 para el siguiente problema y 0 (t) =
y 2 + 2ty t2
0 < t < 1;
y(0.1) =
0.01 . 1.9
1.2 1
yi
0.8 0.6 0.4 0.2 0 0
0.2
0.4
0.6
0.8
1
y-Euler y-RK
ti
y exacta
Figura 8: Aproximaci´on por Euler y RK4. La soluci´on exacta es y(t) = valores se han unido por rectas.
t2 2−t
y los
Otra alternativa ha sido trabajar con el problema integral que se obtiene al integrar respecto a t ambos lados de la igualdad de la ecuaci´on (3), lo que da lugar a Z t y(t) = y(a) + f (s, y(s))ds. a
Al aproximar num´ericamente la integral por medio de las f´ormula de integraci´on de Newton-Cotes se obtienen distintos m´etodos conocidos con el nombre de m´etodos de multipaso. Probar la convergencia tanto de los m´etodos RK, como los de multipaso no es sencillo. En 1956, Dahlquist demostr´o que un m´etodo es convergente si es 20
consistente y estable. Consistente significa que el error local R(t) del m´etodo va a cero cuando h tiende a cero y es estable si el error global permanece acotado y es del mismo orden que el error local. Actualmente la investigaci´on sigue abierta en el caso de las ecuaciones diferenciales algebraicas que aparecen en las aplicaciones de la teor´ıa de control y del c´alculo de variaciones, ver [2], y en problemas de sistemas din´amicos. En M´exico problemas num´ericos en EDO los han trabajado C. Vargas del Cinvestav, F. Solis del CIMAT, R. England, J. L´opez E. y L. Velasco de la FC, entre otros.
6.1
Ecuaciones en derivadas parciales
Los problemas de ecuaciones en derivadas parciales son el fundamento matem´atico de muchos modelos que se usan en la f´ısica, ingenier´ıa, biolog´ıa y recientemente en la econom´ıa y las finanzas. En particular, la mec´anica del medio continuo en dos y tres dimensiones se formula matem´aticamente por medio de ecuaciones en derivadas parciales (EDP). El estudio de las ecuaciones lineales data del siglo XIX y se puede encontrar una soluci´on anal´ıtica para casi todas ellas, cuando est´an definidas en regiones sencillas del plano y del espacio como rect´angulos, cubos, c´ırculos, esferas o cilindros. Cuando la regi´on es m´as complicada o la ecuaci´on es no lineal, s´olo queda aproximar la soluci´on num´ericamente. La soluci´on num´erica de las EDP puede obtenerse de muchas formas. Una clase de m´etodos se obtiene al discretizar tanto la regi´on como las ecuaciones en derivadas parciales para obtener un sistema de ecuaciones en diferencias finitas. Estos m´etodos son una generalizaci´on de los m´etodos para EDO y, como se mencion´o en la parte hist´orica, se usan desde 1910. A principio de los a˜ nos sesenta se contaba con una teor´ıa completa para este tipo de m´etodos para las ecuaciones m´as importantes de la f´ısica matem´atica. El m´etodo de Crank-Nicholson es el m´as popular para problemas par´abolicos y el de upwind y Lax-Wendroff para sistemas hiperb´olicos lineales. El principal inconveniente de los m´etodos de diferencias finitas es que no se adaptan bien a problemas definidos en regiones con fronteras curvas. En los a˜ nos cincuenta, apareci´o en la literatura de ingenier´ıa civil y mec´anica la utilizaci´on del m´etodo de elemento finito. Este m´etodo tienen sus antecedentes en los trabajos de Alexander Hrennikoff y de Richard Courant a inicios de 1940. Su fundamento matem´atico es complemente distinto al de los m´etodos en diferencias finitas y se basa en varios resultados del an´alisis funcional y de las EDP desarrolladas previamente por Rayleigh, Ritz y Galerkin. Este procedimiento ha sido muy exitoso en la aproximaci´on num´erica de la soluci´on de ecuaciones el´ıpticas y parab´olicas en dos y tres dimensiones, en particular en la soluci´on de problemas de mec´anica de fluidos o de convecci´on difusi´on. En la Figura 9 se muestra los resultados num´ericos del estudio de la sedimentaci´on de part´ıculas r´ıgidas en un fluido incompresible realizado con el m´etodo de elemento 21
finito; las flechas indican la velocidad del fluido y los colores la presi´on del mismo. En el caso de los sistemas hiperb´olicos como la ecuaci´on de onda, se usan m´etodos en diferencias finitas especialmente dise˜ nados para adaptarse a soluciones con saltos o con ondas de choque como en los problemas de flujo supers´onico. Para conocer m´as sobre este tema, ver [14] y [4]. A continuaci´on se mencionan algunos grupos y personas de los muchos que han trabajado en este campo en M´exico. Entre los investigadores de la UNAM se cuentan a G. Alduncin y al grupo de I. Herrera en el Instituto de Geof´ısica; en Ciencias de la Atm´osfera el grupo de modelos clim´aticos de J. Adem y el grupo de modelaci´on atmosf´erica de Y. Skiba; en el Instituto de Ingenier´ıa F. J. S´anchez Sesma, G. Ayala y M. Ch´avez, entre muchos otros. En el IIMAS trabaj´o en m´etodos nodales J.P. Hennart, F. Sabina en mec´anica de s´olidos, A. Minzoni en ondas y tambi´en hay que mencionar las contribuciones de B. Chen mientras trabaj´o en M´exico; en DGSCA P. Gonz´ales Casanova estudia los m´etodos radiales, mientras que P. Barrera de la FC se ha dedicado al estudio de la generaci´on autom´atica de mallas en dos dimensiones mientras que L. Esteva ha trabajado en aplicaciones a biomatem´aticas. En el Departamento de Matem´aticas de la UAM-Iztapalapa E. Baez, H. Ju´arez, A. Nicol´as, P. Saavedra y F. S´anchez han aplicado diversos m´etodos en diferencias finitas y elemento finito a problemas de flujo en medios porosos, mec´anica de fluidos, contaminaci´on atmosf´erica, finanzas y tr´afico vehicular. En el Cinvestav C. Vargas ha hecho contribuciones en aplicaciones a biomatem´aticas; en el Cimat el grupo de J.L. Moreles han trabajado en varios problemas de la industria el´ectrica como tambi´en en dispersi´on de ondas; en la BUAP el grupo de A. Fraguela se dedica a la soluci´on te´orica y num´erica de problemas inversos; F. Castillo del ITESM-Estado de M´exico se dedica a problemas de aplicaci´on a ingenier´ıa y en el IMTA el grupo de A. Aldama trabaja en m´etodos num´ericos aplicados a la hidrolog´ıa.
7
Optimizaci´ on
Optimizar recursos, minimizar costos y maximizar beneficios son cuestiones que siempre est´an presentes en el desarrollo de nuevas tecnolog´ıas y en las ciencias econ´omicoadministrativas. Se llaman problemas de optimizaci´on aquellos en los que se busca determinar una soluci´on ´optima. Los problemas reales involucran miles de variables, por lo que dif´ıcilmente se pueden resolver sin el uso de una computadora. La clasificaci´on de los problemas de optimizaci´on puede ser muy amplia aunque nos limitaremos al caso en que la funci´on F a optimizar, llamada funci´on objetivo, sea lo suficientemente regular para usar el c´alculo diferencial. Si la soluci´on se busca en todo Rn se dice que el problema es sin restricciones; en este caso el procedimiento a seguir es directo. Se buscan primero los puntos en los que las derivadas parciales se hacen 22
cero, llamados puntos cr´ıticos, y posteriormente se comprueba si son un m´aximo o un m´ınimo. Esto implica resolver un sistema no lineal, de n ecuaciones con n inc´ognitas. En la mayor´ıa de los casos es un problema dif´ıcil ya que es poco probable que se pueda obtener una soluci´on expl´ıcita. La u ´nica alternativa es aproximar la soluci´on en forma num´erica a trav´es de m´etodos iterativos como el m´etodo de Newton. En la soluci´on num´erica se prefiere trabajar directamente con el problema de optimizaci´on en lugar de resolver el sistema no lineal ya que este enfoque permite −−→ usar la informaci´on del gradiente de la funci´on: ∇F , ver glosario, para acercarse m´as r´apido a la soluci´on. Todos los m´etodos son de tipo iterativo, el m´as sencillo es el de m´aximo descenso; la idea detr´as de este procedimiento es muy geom´etrica: sup´ongase → que se busca el m´ınimo de una funci´on y se toma como punto inicial a − x0 , como → − −−→ d0 = −∇F (x0 ) es la direcci´on en la que disminuye m´as el valor de la funci´on F , se − → → traza la recta que pasa por − x0 con direcci´on d0 y se determina el punto, a lo largo − → → de esta recta, en el que F toma el valor m´as peque˜ no. Este punto x1 = − x0 + t∗ d0 se escoge como el siguiente punto de la iteraci´on y se repite el proceso hasta que la → norma del gradiente en el punto − xi sea menor que una tolerancia dada, lo que implica, − → en principio, que xi es una buena aproximaci´on a la soluci´on. Para seleccionar un m´etodo iterativo se compara costo, n´ umero de operaciones, y rapidez de convergencia. Se dice que un m´etodo iterativo converge con orden p, si el error absoluto en la i-´esima iteraci´on satisface que ei ≤ Kepi−1 , con K > 0. El m´etodo de descenso pronunciado es de orden uno, por lo que el m´etodo es lento pero barato porque s´olo requiere de la evaluaci´on del gradiente en cada iteraci´on. Una alternativa m´as r´apida, con convergencia cuadr´atica, es el m´etodo de Newton que es una adaptaci´on del m´etodo de Newton-Raphson para resolver ecuaciones no lineales. Para problemas cuadr´aticos, Newton converge en una iteraci´on. El procedimiento es similar al de descenso pronunciado con la diferencia que el vector direcci´on se obtiene al resolver en cada iteraci´on un sistema de ecuaciones lineales, cuya matriz asociada es la matriz Hessiana de F, ver glosario. Como puede observarse el costo de cada iteraci´on de Newton es alto, pues hay que evaluar tanto el gradiente como la matriz Hessiana en cada iteraci´on. Para conocer m´as sobre optimizaci´on sin restricciones se sugiere la lectura de [15], uno de los autores es mexicano y se ha distinguido por sus contribuciones a este campo. Hay que mencionar que en los u ´ltimos a˜ nos han aparecido otras t´ecnicas de tipo heur´ıstico como los algoritmos gen´eticos que compiten en algunos casos con los m´etodos aqu´ı mencionados. Un m´etodo es heur´ıstico cuando no se puede garantizar de antemano la convergencia al ´optimo, pero cuyo desempe˜ no ha dado prueba de su efectividad.
23
Optimizaci´ on con restricciones Otro tipo de problemas de optimizaci´on son aquellos en los que se da un conjunto Ω y se desea determinar los puntos de Ω en los que F alcanza su valor m´aximo o m´ınimo. Por lo general Ω est´a definido por un conjunto de restricciones expresadas en forma de ecuaciones: gi (x) = 0 o de desigualdades gi (x) ≤ 0 para i = 1, . . . , m. A estos problemas se les conoce como problemas de optimizaci´on con restricciones. Si F y gi son lineales se tiene un problema de programaci´on lineal. Desde 1947, se conoce un m´etodo iterativo para resolverlo que es el m´etodo simplex creado por Dantzig. Este m´etodo fue el campe´on por muchos a˜ nos hasta que Karmarkar, en 1984, creara un m´etodo muy distinto; sin embargo hasta ahora no se ha podido demostrar la superioridad de uno sobre el otro. Cuando tanto la funci´on objetivo F como las restricciones gi son no lineales y ´estas son ecuaciones se puede aplicar el m´etodo de multiplicadores de Lagrange para obtener un sistema de ecuaciones no lineales de n + m variables. En caso de que las gi sean desigualdades, no se cuenta con un m´etodo que sea la panacea; se tiene, en cambio, m´etodos ad-hoc que aprovechan las caracter´ısticas del conjunto Ω, y las de la funci´on objetivo. Entre los m´etodos de car´acter general m´as populares est´an el de gradiente proyectado, gradiente reducido generalizado, los m´etodos de penalizaci´on y el m´etodo cuadr´atico secuencial. Para saber m´as, consultar [13]. Hay otros problemas de optimizaci´on igualmente importantes, pero que, por falta de espacio, no se incluyen. Problemas de c´alculo de variaciones en los que la funci´on objetivo es un funcional; por ejemplo, una integral y el m´ınimo es una funci´on que pertenece a un espacio de dimensi´on infinita. Otros problemas aparecen en teor´ıa de control. Para saber m´as, consultar [17]. En M´exico han trabajado problemas de optimizaci´on: I. Gitler del Cinvestav, G. Calvillo, Z. Parada del ITAM y C. Gigola e I. Romero de la ESFM del I.P.N.; J. Goddard, M. A. Gutierrez, S. de Los Cobos y B. P´erez en la UAM-Iztapalapa y C. Barr´on en la UAM-Cuajimalpa. En la UNAM, los grupos de P. Barrera y J. L´opez E. de la FC, D. Romero del IM, A. Velasco Levy, S. G´omez y Lu´ıs Morales del IIMAS; R. R´ıos de la UANL y P. Flores de UNISON, entre otros.
8
Perspectivas a futuro
Es dif´ıcil prever de antemano el futuro de un campo cuyo avance depende tanto del desarrollo tecnol´ogico; nuevos problemas aparecer´an conforme se creen nuevas tecnolog´ıas. En 1950, nadie pod´ıa imaginar el papel que desempe˜ nan, en la vida cotidiana de hoy en d´ıa, las computadoras y los instrumentos electr´onicos. Sin embargo, es de esperarse que el cambio clim´atico, la astronom´ıa, que contin´ ua siendo
24
una fuente inagotable de problemas abiertos y el desarrollo de nuevos materiales, entre otros muchos temas, demanden de m´etodos num´ericos m´as r´apidos y m´as precisos que aprovechen el creciente incremento de capacidad computacional. Asimismo, se espera que la modelaci´on matem´atica y la matem´atica num´erica se vuelvan herramientas indispensables en la biolog´ıa, la farmacolog´ıa y la medicina y que ayuden, entre otras cosas, al desarrollo de nuevos medicamentos y hacer una realidad el diagn´ostico m´edico computarizado. Uno de los aspectos que est´a cambiando en la matem´atica num´erica es que cada vez es m´as frecuente que para resolver un mismo problema se usen t´ecnicas h´ıbridas que combinan algoritmos deterministas con estoc´asticos; hasta hace poco, los investigadores se especializaban en s´olo uno de estos enfoques, pero ahora en varios campos como las finanzas, el estudio de la turbulencia, del tr´afico vehicular o la din´amica de poblaciones se usan ambas metodolog´ıas. Esta tendencia se generalizar´a en otros temas, lo que exigir´a una formaci´on m´as amplia y flexible por parte del matem´atico num´erico. Otros aspectos importantes en los que se har´a m´as ´enfasis en el futuro en el dise˜ no de algoritmos son la paralelizaci´on y la adaptabilidad. Lo primero implica que a nivel de programaci´on se explote, al cien por ciento, la posibilidad de realizar distintas tareas al mismo tiempo con distintos procesadores; la adaptabilidad implica que un algoritmo sea capaz de ajustar el valor de sus par´ametros dependiendo de lo que mejor convenga al tipo de problema que se est´a resolviendo. Cada vez m´as, los paquetes computacionales que tienen ´exito son aquellos que tienen mayor aplicabilidad y capacidad de adaptarse autom´aticamente a los requerimientos de cada problema con el fin de aproximar la soluci´on de la manera m´as precisa en el menor tiempo posible. Por u ´ltimo, es de esperarse, que a medida que sigan evolucionando los instrumentos de c´alculo y el grado de complejidad de los problemas a resolver, seguir´an ampli´andose los campos de acci´on de la matem´atica num´erica.
9
Glosario 1. N´ umero de condici´on de una matriz cuadrada invertible se define por: cond(A) = ||A||2 ||A−1 ||2 , con ||A||2 igual a la ra´ız cuadrada del valor propio m´as grande en valor absoluto de la matriz At A. 2. Diagonal de una matriz: los elementos Aii , con i = 1, . . . , n, de una matriz cuadrada. 3. Funci´on Lipschitz. Una funci´on f es Lipschitz en el conjunto [a, b] si existe un real L > 0 tal que para cualquier x, y ∈ [a, b] |f (x) − f (y)| ≤ L|x − y|. 25
4. ∇F (x0 ) = ( ∂F∂x(x10 ) , . . . , ∂F∂x(xn0 ) ). 5. Inversa de una matriz: Dada una matriz cuadrada A, la matriz inversa de A se denota por A−1 y es la matriz que satisface: AA−1 = A−1 A = I con I la matriz identidad. 6. Matriz Hessiana H : Matriz cuyo elemento Hij =
∂2F . ∂xi ∂xj
7. Matriz identidad: es una matriz cuadrada con unos en la diagonal y cero en el resto de los elementos. 8. Matriz sim´etrica es una matriz que es igual a su transpuesta. 9. Matriz positiva definida es aquella que satisface xt Ax > 0 para todo vector x 6= 0. 10. Matriz triangular superior: una matriz cuadrada cuyos elementos abajo de la diagonal son cero. Matriz triangular inferior: una matriz cuadrada cuyos elementos arriba de la diagonal son cero. 11. Rn : es el conjunto de vectores con n componentes en los reales. 12. La matriz transpuesta de una matriz A se denota por At y Atij = Aji , para i, j = 1, . . . , n.
References [1] P. Barrera, V. Hern´andez y C. Dur´an. El ABC de los splines. Aportaciones Matem´aticas. Textos No. 9. Publicaciones SMM. 1996. [2] K.E. Brenan, S.L. Campbell y L.R. Petzold. Numerical Solution of Initial Value Problems in Differential Algebraic Equations. Classics in Applied Mathematics. 14. SIAM. 1996. [3] J.P. Collette. Historia de las matem´aticas I y II. Tercera edici´on. Siglo XXI Editores. 1998. [4] P. G. Ciarlet. The Finite Element Method for Elliptic Problems. North-Holland. 1978. [5] J. Demmel. Applied Numerical Linear Algebra. SIAM. 1997. [6] G.E. Forsythe. Pitfalls in computation, or why a math book isn’t enough. The American Monthly, Vol. 77, No. 9. 931-956. 1970. 26
[7] G. Golub y C. Van Loan. Matrix Computations. Tercera edici´on. John Hopkins University Press. 1996. [8] H. Goldstein. A History of Numerical Analysis: From 16th through 19th Century. Springer Verlag. 1977. [9] H. Goldstein y J. V. Neumann. On the principles of large scale computing machines. Proceedings of the Summer Research Institute on the Legacy of John Von Neumann. Vol 50, AMS. 1990. [10] A. Greenbaum. Iterative Methods for Solving Linear Systems. SIAM. 1997. [11] D.B. Hern´andez. An´alisis Num´erico. III Coloquio del Departamento de Matem´aticas del Cinvestav. 1983. [12] H. Madrid. El ´algebra lineal num´erica del PageRank de Goggle. Memorias de la Sociedad Matem´atica de M´exico 2006. Aportaciones Matem´aticas. Comunicaciones 36. 33-52. 2006. [13] J.M. Martinez y S. A. Santos. Metodos Computacionais de Optimizacao. Notas del curso impartido en el 20 Coloquio Brasileiro de Matem´aticas. IMPA. 1995. [14] K. W. Morton y D.F. Mayers. Numerical Solution of Partial Differential Equations. Second Edition. Cambridge University Press. 2005. [15] Nocedal, J. y S.W. Wright. Numerical Optimization. Springer. 1999. [16] M. Overton. C´omputo num´erico con aritm´etica en punto flotante IEEE. Aportaciones Matem´aticas. Textos Nivel Medio. Num. 19. SIAM y SMM. 2002. [17] D.A. Pierre. Optimization Theory with Applications. Dover. Second Edition. 1986. [18] C. Vargas. Soluci´on num´erica de ecuaciones diferenciales. IV Coloquio del Departamento de Matem´aticas del Cinvestav. 1985. [19] E. Trabulse. Historia de la ciencia en M´exico. Fondo de Cultura Econ´omica. 1997. [20] L.N. Trefethen. Numerical Analysis. 2006 ( por aparecer en el Princeton Companion to Mathematics. Editado por T. Gowers y J. Barrow-Green y publicado por Princeton University Press).
27
t=0
t=1
t=2
t=3
Figura 9: Sedimentaci´on de 100 part´ıculas r´ıgidas en un fluido viscoso incompresible. Cortes´ıa de H. Ju´arez
28