Introducci´on a la F´ısica Computacional: Problemas de F´ısica resueltos num´ericamente Facultad de Ciencias, Universidad de Oviedo Curso 2014-15 Julio M. Fern´andez D´ıaz, Rosario D´ıaz Crespo

Introducci´on Estos apuntes de aplicaciones en F´ısica resueltos num´ericamente abordan diversos problemas de F´ısica Computacional, con un doble objetivo: Introducir ejemplos concretos de problemas f´ısicos que deben resolverse, y que tienen relaci´on directa con problemas matem´aticos determinados (ceros de funciones, minimizaci´on, etc.). Esto implica realmente introducir familias de problemas similares que los estudiantes se encontrar´an a lo largo de sus estudios. Desarrollar m´etodos num´ericos sencillos, que pueden entenderse f´acilmente a nivel de primero de grado, que sirven como introducci´on a los m´etodos m´as avanzados que los estudiantes abordar´an en cursos superiores. Estos m´etodos sencillos se comparan (en cuanto a resultados) con los que ofrece el m´odulo scipy, de Python, que son evidentemente los recomendados de manera general. Los seis problemas que se resolver´an son los siguientes: 1. Distribuci´on de Boltzmann para las energ´ıas de las mol´eculas de un gas: M´etodo de Montecarlo. 2. Ecuaci´on de estado de un gas: Ra´ıces de ecuaciones no lineales. 3. Determinaci´on num´erica del per´ıodo de un p´endulo de amplitud cualquiera: Cuadraturas (integrales) en una dimensi´on. 4. Movimiento unidimensional: ca´ıda de un paracaidista: Problema diferencial de valores iniciales. 5. C´alculo num´erico de centros de masa y de momentos de inercia: Cuadraturas (integrales) en dos dimensiones. 6. Velocidad de propagaci´on de ondas superficiales en un l´ıquido: C´alculo de derivadas y miminimizaci´on de funciones. Obviamente existen muchos m´as problemas de F´ısica y m´etodos num´ericos de inter´es, pero la falta de tiempo y el nivel de los estudiantes impiden una extensi´on mayor de estos apuntes.

Problema 1 Distribuci´on de Boltzmann para las energ´ıas de las mol´eculas de un gas (M´etodo de Montecarlo) El m´etodo de Montecarlo es quiz´a el algoritmo m´as fruct´ıfero de la F´ısica Computacional. En sus diferentes variantes permiti´o la simulaci´on de muchos fen´omenos f´ısicos que tienen un componente aleatorio: absorci´on y transmisi´on de part´ıculas gamma en una laja de material, movimiento de part´ıculas en un gas o l´ıquido, agregaci´on de part´ıculas, etc. La base del m´etodo es la obtenci´on de n´umeros aleatorios que cumplen cierta distribuci´on. Un ordenador no puede por principio generar n´umeros exactamente aleatorios puesto que es una m´aquina determinista. Sin embargo se dise˜nan m´etodos para obtener n´umeros que cumplan ciertas propiedades de aleatoriedad que se denominan entonces pseudoaleatorios. El m´etodo de Montecarlo puede combinarse con cualesquiera otros m´etodos de la F´ısica Computacional y del C´alculo Num´erico. Tiene por tanto una gran importancia.

1.1. Objetivos Analizar la evoluci´on hacia el equilibrio de un sistema termodin´amico simple, formado por mol´eculas de un gas que sufren choques en los que intercambian su energ´ıa. En la resoluci´on del presente problema se utiliza el m´etodo de Montecarlo, de tanta importancia en la simulaci´on de problemas de todo tipo en la actualidad.

1.2. Conocimientos previos Estad´ıstica: histogramas.

1

Universidad de Oviedo

Departamento de F´ısica

Termodin´amica: concepto de equilibrio. F´ısica estad´ıstica: distribuci´on de Boltzmann.

1.3. Planteamiento del problema El comportamiento termodin´amico de un sistema es debido fundamentalmente a que est´a formado por numerosos componentes (en un gas, las mol´eculas del mismo). En general habr´a un n´umero muy grande de part´ıculas, del orden de 1023 o m´as en un sistema macrosc´opico. Los sistemas inicialmente en desequilibrio, si se dejan evolucionar, tienden espont´aneamente al equilibrio termodin´amico. La evoluci´on termodin´amica es siempre hacia estados m´as desordenados que los originales, esto es, de mayor entrop´ıa. Un sistema macrosc´opico como un volumen dado de un gas contiene numerosas mol´eculas en continuo movimiento. Cada una de las mol´eculas presentes posee cierta energ´ıa cin´etica (y de otros tipos) en cada momento, pero a lo largo del tiempo, debido a los choques intermoleculares las part´ıculas van intercambiando su energ´ıa. Sin embargo, debido a la existencia de la colectividad, el n´umero de part´ıculas con una energ´ıa dada no puede ser cualquiera.

1.4. Distribuci´on te´orica de energ´ıas en un conjunto de part´ıculas en interacci´on La temperatura del gas es directamente proporcional a la energ´ıa cin´etica media de las mol´eculas que lo componen. Esto significa que si la temperatura es alta las mol´eculas tienen en media grandes energ´ıas. Sin embargo siempre pueden existir mol´eculas con energ´ıas tan grandes como queramos. Existe por tanto una distribuci´on estad´ıstica (pr´acticamente continua debido al gran n´umero de mol´eculas) de velocidades y de energ´ıas moleculares. Para la energ´ıa esa distribuci´on se denomina de Boltzmann. La misma se obtiene dejando evolucionar al sistema, a lo largo del tiempo, y los choques intermoleculares se encargan de que se alcance esa distribuci´on (como l´ımite cuando t → ∞), siempre, independientemente de la distribuci´on inicial de energ´ıas. Te´oricamente, se puede demostrar que si la energ´ıa media es Em la funci´on de densidad correspondiente es:   dn N E , = exp − dE Em Em

donde N es el n´umero total de part´ıculas presentes. A partir de esta expresi´on se puede demostrar f´acilmente que existir´a en todo momento una fracci´on   ∆n E 1 exp − ∆E ≈ N Em Em de mol´eculas con energ´ıas entre E − 12 ∆E y E + 12 ∆E. Las mol´eculas no ser´an siempre las mismas, debido a los continuos choques. Nosotros obtendremos esa distribuci´on de energ´ıas utilizando un m´etodo num´erico, no una demostraci´on te´orica. Introducci´on a la F´ısica Computacional

2

Universidad de Oviedo

Departamento de F´ısica

1.5. El m´etodo de Montecarlo de simulaci´on num´erica Este m´etodo consiste en simular num´ericamente fen´omenos f´ısicos (y de otras ciencias) usando variables aleatorias cuya distribuci´on se ajusta al fen´omeno tratado. Fue inventado por Von Neumann y Ulam para c´alculos dif´ıciles que deb´ıan realizar en el dise˜no de artefactos nucleares, y bautizado as´ı por su analog´ıa con los juegos de azar de los casinos. Antes de poder aplicar el m´etodo de Montecarlo, es interesante recordar las nociones de probabilidad y de variable aleatoria. Desde un punto de vista estrictamente matem´atico tales cuestiones son dif´ıciles y no podemos entrar a fondo en el tema, conform´andonos con aproximaciones que se entiendan bien.

1.6. Probabilidad y variables aleatorias Un experimento de lanzamiento de monedas sirve para explicar el concepto de probabilidad: n´umero de casos favorables dividido por el de casos posibles, una vez que tiramos la moneda un n´umero infinito de veces (poblaci´on de tiradas). Aqu´ı llegamos a que la probabilidad es un n´umero entre 0 y 1 que asignamos a un suceso. Adem´as, todo el conjunto de sucesos posible tiene una probabilidad 1 de ocurrir (suceso seguro). Puesto que en la pr´actica es imposible realizar infinitas tiradas, se realiza un n´umero finito de ellas, experiencia denominada muestra. Se define entonces variable aleatoria, como una magnitud que puede tomar todos los valores iguales a los sucesos y s´olo esos, y adem´as la probabilidad de que lo haga es igual a la probabilidad de ocurrencia del suceso. En nuestro ejemplo de la moneda, una variable aleatoria puede tomar dos valores discretos: cara o cruz. En el caso planteado se tiene una variable aleatoria discreta. Si el conjunto de valores posibles es continuo la variable aleatoria se denomina continua.

´ 1.6.1. Numeros aleatorios Un n´umero aleatorio es un n´umero obtenido por azar, pero no de cualquier manera sino siguiendo una ley de probabilidades (ley de densidad). Un caso f´acil de entender es el de un dado de 6 caras: se pueden obtener los n´umeros 1, 2, 3, 4, 5 y 6. Se tiene entonces una variable aleatoria discreta, que puede tomar s´olo esos valores.

Introducci´on a la F´ısica Computacional

3

Universidad de Oviedo

Departamento de F´ısica

1.6.2. Variables aleatorias discretas Tirar un dado es hacer una realizaci´on de la variable aleatoria: sale uno de esos n´umeros. Si el dado no est´a cargado existe igual probabilidad de que salga cualquier n´umero de 1 a 6, que ser´a posiblemente diferente en cada tirada debido al azar. f es la ley de densidad.

f 1 6

1

2

3

4

5

6

x

1.6.3. Variables aleatorias continuas Se puede extender la definici´on a variables continuas, con distribuciones de probabilidad continuas. En ese caso se trabaja con a´ reas entre dos puntos del eje Ox: R P = ab f (x)dx es la probabilidad de que un n´umero tomado al azar seg´un la distribuci´on f caiga entre a y b.

f

P a

b

x

Dos distribuciones continuas son de mucha importancia: la distribuci´on uniforme y la distribuci´on normal o gaussiana. La distribuci´on uniforme En la pr´actica existe una distribuci´on b´asica: la distribuci´on uniforme entre 0 y 1, etiquetada U [0, 1), que tiene por funci´on de densidad: f

f (x) =



1 0

0≤x>> from numpy import random >>> random.seed(200); print(random.rand(2))

[ 0.94763226 0.22654742] >>> random.seed(200); print(random.rand(2)) [ 0.94763226 0.22654742]

Si no se especifica la semilla se usa t´ıpicamente el reloj del ordenador (o /dev/urandom si existe en el ordenador). Arrays aleatorios: funciones simples La funci´on random.rand(shape) devuelve un array con la forma dada relleno de n´umeros U (0, 1). Ejemplo 1.6.2. Generaci´on de una matriz (2,3) de n´umeros aleatorios con distribuci´on uniforme entre 0 y 1: >>> from numpy import random >>> print(random.rand(2,3)) [[ 0.95344434, 0.76328461, 0.26307595], [ 0.606801 , 0.16873616, 0.06620183]] random.randn(shape) devuelve un array con la forma dada relleno de n´umeros con distribuci´on N(0, 1). Ejemplo 1.6.3. Generaci´on de una matriz (3,2) de n´umeros aleatorios con una distribuci´on normal con media 0 y desviaci´on t´ıpica 1: >>> print(random.randn(3,2)) [[-1.01948701, 0.27342484], [-0.02932765, 1.06019774], [ 0.85948562, 2.03222582]] random.normal(m, s, shape) devuelve un array con la forma dada relleno de n´umeros N(m, s). Ejemplo 1.6.4. Generaci´on de un vector de 2 elementos de n´umeros aleatorios con una distribuci´on normal de media 10 y desviaci´on estandar 0.1: >>> print(random.normal(10,.1,2,)) [ 9.90764775, 9.97226389] Introducci´on a la F´ısica Computacional

6

Universidad de Oviedo

Departamento de F´ısica

random.uniform(a, b, shape) devuelve un array con la forma dada relleno de n´umeros U (a, b). Ejemplo 1.6.5. Generaci´on de un vector de 2 elementos de n´umeros aleatorios con distribuci´on uniforme entre 10 y 20: >>> random.uniform(10,20,2,) array([ 16.72835576, 10.32085779]) >>> print(random.uniform(10,20,(2,2))) [[ 14.73801182, 18.45998782], [ 14.64104923, 11.4765427 ]] Arrays aleatorios: distribuciones Se obtienen n´umeros aleatorios que siguen otras distribuciones usando un cambio de variable aleatoria. El tema es demasiado espec´ıfico para analizarlo aqu´ı. El m´etodo, a grandes rasgos, consiste en generar n´umeros aleatorios distribuidos uniformemente en [0, 1) y mediante una transformaci´on obtener el n´umero aleatorio que nos dar´ıa otra distribuci´on. En numpy.random se incluyen algunas funciones que nos proporcionan el n´umero buscado seg´un distribuciones importantes para los estad´ısticos: random.binomial(n, p) binomial. random.chisquare(df) chi cuadrado. random.poisson(lam) poisson. random.lognormal(mean, sigma) lognormal. random.standard t(df) t de Student. random.f(dfnum, dfden) F de Fisher-Snedecor.

1.6.6. Ejercicios sobre generaci´on de aleatorios Ejercicio 1.6.1. Prepare un programa que simule un sorteo del juego de Euromillones que consiste en sacar una serie (denominada ‘n´umeros’) de 5 n´umeros de un total de 50 y de otra (denominada ‘estrellas’) de 2 n´umeros de un total de 11. Ejercicio 1.6.2. Prepare un programa que simule la tirada de dos dados, obteniendo la suma del resultado de cada dado como valor aleatorio.

Introducci´on a la F´ısica Computacional

7

Universidad de Oviedo

Departamento de F´ısica

1.7. Desarrollo de un programa que nos analice la evoluci´on energ´etica de un conjunto de part´ıculas Vamos a desarrollar un programa que simule la evoluci´on de un sistema con una distribuci´on cualquiera a partir de los choques equiprobables entre las part´ıculas del sistema, las cuales intercambian tambi´en aleatoriamente su energ´ıa. Se va a trabajar con valores enteros de la energ´ıa. Debido a que la energ´ıa s´olo puede almacenarse en paquetes seg´un la Mec´anica cu´antica, tal suposici´on es bastante realista y nos ahorra diversas complicaciones en la programaci´on y representaci´on de resultados. El uso de energ´ıas enteras tiene un efecto colateral: la distribuci´on de Boltzmann se ajusta peor para valores peque˜nos de la energ´ıa, como podremos observar en nuestras simulaciones. El programa no necesita ninguna funci´on definida por el usuario. En la parte inicial se define el n´umero total de part´ıculas, np, la energ´ıa inicial de todas las part´ıculas, ei (que coincide con la energ´ıa media Em ), y el n´umero de interacciones a efectuar, ni. Para efectuar el dibujo se necesita tambi´en una energ´ıa m´axima para las abscisas, emax, que hacemos igual a 10*ei, y cada cuantas interaciones (nis) haremos el gr´afico. La energ´ıa de cada part´ıcula se guarda en un vector entero, ep, inicializado a ei. Asimismo, se usa otro vector entero, ne, con el n´umero de part´ıculas que tienen una energ´ıa dada. Para realizar las interacciones se debe primero determinar qu´e dos part´ıculas del sistema interaccionan. Puesto que la probabilidad de interacci´on es la misma para todas las part´ıculas podemos usar la funci´on randint(k) del m´odulo random de numpy, que nos devuelve un valor aleatorio entero entre 0 y k-1 con igual probabilidad. Una vez obtenidas aleatoriamente de esa manera dos part´ıculas (diferentes), se suman sus energ´ıas, y se vuelven a repartir tambi´en aleatoriamente. Tambi´en hay que actualizar el vector ne en cada interacci´on. Este proceso se realiza ni veces. La salida de resultados se realiza en un bloque de c´odigo con gr´aficos de matplotlib, cada nis interacciones. No es necesario vectorizar el c´odigo demasiado, puesto que la salida gr´afica es lenta y no se gana mucho en eficiencia vectorizando el c´odigo. # -*- coding: utf-8 -*# Distribuci´ on de energ´ ıa de Boltzmann from numpy import * from numpy.random import randint import matplotlib.pyplot as plt np = 100000 ei = 5 ni = 800000 nis = 10000

# # # #

n´ umero total de part´ ıculas energ´ ıa inicial de cada part´ ıcula n´ umero de interacciones n´ umero de interacciones entre salida de resultados

Introducci´on a la F´ısica Computacional

8

Universidad de Oviedo

emax = 10*ei

Departamento de F´ısica

# energ´ ıa m´ axima para dibujar

ep = ei*ones(np, int) # se inicializa la energ´ ıa de cada part´ ıcula ne = zeros(emax+1) # contador de part´ ıculas para cada energ´ ıa ne[ei] = np # todas las part´ ıculas tienen energ´ ıa inicial ’ei’ # distribuci´ on te´ orica e = arange(emax+1) nt = 1.0/ei*exp(-e/float(ei))

# la energ´ ıa media es ’ei’

for i in range(ni): # salida de resultados cuando ’i’ sea multiplo de ’nis’ if i % nis == 0: plt.clf() nef = ne/float(np) plt.bar(e, nef, width=0.6, color=’b’) plt.plot(e, nt, ’-r’) fmax = 2*max(nt[0], max(nef)) plt.axis([0, emax, 0, fmax if fmax < 1.0 else 1.0]) plt.xlabel(u’Energ´ ıa’) plt.ylabel(u’Fracci´ on de part´ ıculas con energ´ ıa dada’) plt.title(u’interacci´ on n´ umero ’+str(i)) plt.draw() plt.show(block=False) # se sortean dos part´ ıculas diferentes ip1 = randint(np) while True: ip2 = randint(np) if ip1 != ip2: break et = ep[ip1] + ep[ip2] # energ´ ıa total de ambas particulas if et == 0: continue ne[ep[ip1]] -= 1 ne[ep[ip2]] -= 1 # la energia ’et’ se reparte entre ambas particulas aleatoriamente # las part´ ıculas con e > emax se asocian a ’emax’ exactamente e1 = randint(et+1) if e1 > emax: e1 = emax e2 = et-e1 if e2 > emax: e2, e1 = emax, et-emax ep[ip1] = ep[ip2] = ne[e1] += ne[e2] +=

e1 e2 1 1

Introducci´on a la F´ısica Computacional

9

Universidad de Oviedo

Departamento de F´ısica

En el progama ejemplo se toma ei = 5. Se deja como ejercicio probar otros casos con diferentes n´umeros totales de part´ıculas y energ´ıas iniciales (por ejemplo, ei = 20).

1.8. Ejemplos adicionales y de ´ındole pr´actica Distribuci´on de velocidades de Maxwell-Boltzmann.

Introducci´on a la F´ısica Computacional

10

Problema 2 Ecuaci´on de estado de un gas (Ra´ıces de ecuaciones) La obtenci´on de ra´ıces de ecuaciones (tambi´en denominada obtenci´on de ceros de funciones) es una parte b´asica de cualquier biblioteca de m´etodos num´ericos. Existen muchos m´etodos, m´as o menos robustos, para obtener ra´ıces de ecuaciones. Algunos de ellos est´an preparados para funciones concretas (por ejemplo polinomios). Se suelen distinguir los problemas de una dimensi´on de los de varias dimensiones, puesto que aquellos son bastante m´as simples de abordar. Existen incluso m´etodos robustos (aunque algo lentos) como la bisecci´on, dif´ıciles de implementar en varias dimensiones. En los casos de varias dimensiones se suelen utilizar m´etodos de gradiente (derivada) en los que la funci´on se ‘lleva’ hacia cero por una pendiente descendiente.

2.1. Objetivos Estudiar la ecuaci´on de estado de un gas no ideal (concretamente la ecuaci´on de Van der Waals), comparando los resultados obtenidos con el caso ideal. Para tratar el problema se necesitar´a resolver ecuaciones no lineales, y se ver´an m´etodos de resoluci´on de este tipo de ecuaciones. Se plantea una soluci´on simple por el m´etodo de bisecci´on y luego se aborda la soluci´on usando la funci´on fsolve de scipy.

2.2. Conocimientos previos Matem´aticas: ecuaciones en el cuerpo de los n´umeros reales. Termodin´amica: ecuaci´on de estado.

11

Universidad de Oviedo

Departamento de F´ısica

Tabla 2.1: Constantes de Van der Waals de varios gases. gas CO2 O2 butano

a

b

(Pa m6 mol−2 )

(10−5 m3 mol−1 )

0.3640 0.1378 1.466

4.267 3.183 12.26

2.3. Ecuaciones de estado de los gases reales Las variables termodin´amicas que describen el estado de un sistema no pueden tomar valores cualesquiera: una vez fijadas algunas, otras obligatoriamente toman ciertos valores definidos a partir de ecuaciones que describen el estado del sistema. Los casos m´as simples de ecuaciones de estado las cumplen los gases. Es sobradamente conocida la ecuaci´on de estado de un gas ideal: PV = nRT, donde P es la presi´on, V el volumen ocupado por el gas, n el n´umero de moles presentes, T la temperatura absoluta, y R una constante de valor 8.314472 J K−1 mol−1 . En ocasiones se trabaja con el volumen molar: v = V /n, de tal manera que tendr´ıamos: Pv = RT.

(2.1)

La ecuaci´on de los gases ideales es v´alida para gases reales cuando las presiones son peque˜nas y/o las temperaturas no son muy bajas. Para gases reales se han propuesto diversas ecuaciones m´as exactas en condiciones no ideales. Una muy conocida es la ecuaci´on de Van der Waals:  a (2.2) P + 2 (v − b) = RT, v donde a y b son dos par´ametros que dependen del gas que estemos analizando. En la tabla 2.1 se exponen valores para diversos gases. Otra de las ecuaciones de estado no demasiado complejas es la de Redlich-Kwong: RT a , (2.3) P= −√ v−b T v(v + b) donde a y b son dos par´ametros que dependen del gas que estemos analizando. En la tabla 2.2 se exponen valores para diversos gases. Las tres ecuaciones planteadas tienen en com´un que la presi´on, P, aparece linealmente en ellas. En la ecuaci´on de los gases ideales y en la de Van der Waals la temperatura, T , tambi´en aparece linealmente, mientras que el volumen molar, v, aparece s´olo linealmente en la ecuaci´on de los gases ideales. Esto significa, por ejemplo, que para obtener el volumen molar que ocupa un gas dado a partir de las otras variables termodin´amicas se necesita resolver una ecuaci´on no lineal. Introducci´on a la F´ısica Computacional

12

Universidad de Oviedo

Departamento de F´ısica

Tabla 2.2: Constantes de la ecuaci´on de Redlich-Kwong de varios gases. gas CO2 O2 butano

a

b

(Pa K0.5 m6 mol−2 )

(10−5 m3 mol−1 )

6.447 1.741 29.016

2.963 2.208 8.068

2.4. M´etodos num´ericos de c´alculo de ceros de funciones Existen muchos m´etodos para resolver, en el cuerpo de los reales, la ecuaci´on f (x) = 0.

(2.4)

Analizaremos el m´etodo de bisecci´on que es seguro (si se cumplen ciertas condiciones b´asicas, siempre se encuentra una ra´ız). Este m´etodo tiene una convergencia lineal, de tal manera que se llega lentamente a una ra´ız con una aproximaci´on dada. Otros m´etodos, como el de Newton-Raphson, que usa tambi´en la derivada de la funci´on de la que se quiere la ra´ız, tienen convergencia m´as r´apida. El m´etodo de bisecci´on, basado en el teorema de Rolle del An´alisis matem´atico, es usualmente lento, pero bajo ciertas condiciones nos asegura encontrar la ra´ız buscada, permiti´endonos tambi´en controlar el error con que se ha encontrado la misma. El teorema de Rolle nos dice que toda funci´on continua en un intervalo con valores opuestos en signo en los extremos del mismo debe poseer al menos un cero en ese intervalo. La aplicaci´on directa de este teorema, con la subdivisi´on sucesiva del intervalo nos lleva a una ra´ız indefectiblemente. El algoritmo es el siguiente: 1. Se parte de dos abscisas x1 y x2 tales que f (x1 ) f (x2 ) < 0. 2. Se calcula x¯ = 12 (x1 + x2 ). 3. Si f (x) ¯ = 0 o´ |x2 − x1 | < ε se toma x¯ como ra´ız y se acaba. 4. Si f (x) ¯ f (x2 ) < 0 se hace x1 = x. ¯ En caso contrario se hace x2 = x. ¯ 5. Se vuelve al punto 2. El valor de ε es el error absoluto que queremos que tenga el valor encontrado para la ra´ız en relaci´on al valor exacto, por ejemplo 10−6 . Si se parte de dos valores iniciales dados x1 = a y x2 = b, el error del m´etodo en la iteraci´on n es menor que 2−n |b − a|. De aqu´ı que si se quiere una soluci´on con un error absoluto dado menor que una tolerancia dada xtol , se debe usar un n´umero de iteraciones dado por: nmin = ⌈log2 (|b − a|/xtol )⌉. Introducci´on a la F´ısica Computacional

13

Universidad de Oviedo

Departamento de F´ısica

En el m´etodo de bisecci´on es fundamental partir de dos abscisas x1 y x2 que cumplan f (x1 ) f (x2 ) < 0. Se necesita que la ra´ız est´e ‘encerrada’ (en ingl´es bracketing), y est´a es la mayor dificultad del m´etodo: buscar dos abscisas adecuadas (lo que depende del problema a resolver). Los m´etodos que est´an integrados en los paquetes de c´alculo num´erico normalmente est´an preparados para que nos obtengan una ra´ız si existe. Pueden estar basados tanto en bisecci´on como en el m´etodo de Newton-Raphson (u otros parecidos) pero con salvaguardas de tipo num´erico. fsolve de scipy es uno de estos m´etodos. En cualquier caso es fundamental no realizar a ciegas la b´usqueda de los ceros de funciones. En primer lugar puede interesarnos un cero de funci´on en una zona determinada del eje de abscisas, por ejemplo porque sea la zona de inter´es f´ısico. No s´olo eso, sino que deberemos partir suficientemente cerca de la soluci´on para que el m´etodo converga en ocasiones. Adem´as interesa siempre calcular el valor de la funci´on en el punto determinado como ra´ız como comprobaci´on. Eso significa que es conveniente dibujar (en una amplio rango de valores de la abscisa) la funci´on de la que se quiere determinar el cero, ampliando si es necesario la zona de inter´es en el dibujo. Esto en nuestro caso se puede hacer f´acilmente con matplotlib. Adem´as es conveniente dibujar el eje de abscisas o un mallado (con grid) para ayudarnos en visualizar los cortes con el eje de abscisas. Ejemplo 2.4.1. Ra´ıces de funciones. Use fsolve para determinar todas las ra´ıces reales del polinomio: f (x) = −2.0 + 6.2x − 4.1x2 + 0.7x3 . Soluci´on # -*- coding: utf-8 -*from numpy import * from scipy.optimize import fsolve import matplotlib.pyplot as plt def f (x): return -2.0+x*(6.2+x*(-4.1+0.7*x)) xp = linspace(-1, 4, 200) fx = f(xp) plt.grid() plt.plot(xp, fx, ’-’) plt.show() x1 = fsolve(f, -2.0); print(x1, f(x1)) x2 = fsolve(f, 2.0); print(x2, f(x2)) x3 = fsolve(f, 4.0); print(x3, f(x3)) Ejercicio 2.4.1. Repita el problema anterior con: g(x) = (((x − 1.97)x + 5.0)x − 2.8)x − 5.03 Introducci´on a la F´ısica Computacional

14

Universidad de Oviedo

Departamento de F´ısica

2.5. Desarrollo de un programa que nos determine el volumen molar de un gas En un laboratorio se quiere determinar el volumen molar (en L/mol) de un gas (CO2 concretamente) en funci´on de la presi´on y temperatura de forma que los recipientes que va a contener el gas tenga un tama˜no adecuado en cada condici´on. Adem´as tambi´en interesa determinar c´omo de buena es la ley de los gases ideales en el estudio del comportamiento del gas. Para ello se quieren comparar los resultados que se obtienen con esa ley y con la ecuaci´on de Van der Waals. Concretando se pretende crear una tabla que tenga temperaturas variables desde 300 K hasta 1000 K (con incrementos de 100 K) y, para cada temperatura, presiones variando logar´ıtmicamente entre 1 atm y 100 atm con 6 intervalos por d´ecada. (En el c´alculo hay que tener presente que 1 atm equivale a 101325 Pa.) Por otro lado, es interesante determinar cu´antas veces se invoca a la funci´on de la que se quiere obtener la ra´ız con objeto de estimar la bondad del m´etodo de c´alculo de ra´ıces. Para ello se utiliza una variable global nfun que se inicializa a cero antes de invocar a la funci´on que calcula la ra´ız y que se va incrementando dentro de la funci´on objetivo. En un c´odigo de producci´on (no de ejemplo como el presente) no debe usarse nfun. La ecuaci´on (2.2) debe ponerse de la forma (2.4), aunque con la variable v en vez de x, como por ejemplo: (Pv2 + a)(v − b) − RT v2 = 0. La funci´on fsolve necesita un valor aproximado inicial de la ra´ız buscada. Eso se solventa usando el valor de videal como valor inicial. Por otro lado, biseccion necesita dos valores de volumen que encierren a la ra´ız buscada. Para ello se utilizan 10−6 videal y 10videal como l´ımites de b´usqueda. # -*- coding: utf-8 -*# soluci´ on de la ecuaci´ on de van der Waals para v from numpy import * from scipy.optimize import fsolve # Para la ecuaci´ on de Van der Waals def fVW (v): global nfun # para contar el n´ umero de veces nfun += 1 # que se invoca esta funci´ on return (P*v**2+a)*(v-b)-R*T*v**2 # M´ etodo de bisecci´ on def biseccion (f, a, b, xtol=1.0e-8): fa, fb = f(a), f(b) if fa == 0.0: return a if fb == 0.0: return b Introducci´on a la F´ısica Computacional

15

Universidad de Oviedo

Departamento de F´ısica

if fa*fb > 0: print("Error en bisecci´ on") exit(1) n = 1+int(log(abs(b-a)/xtol)/log(2.0)) for i in range(n): x = (a+b)/2.0 fx = f(x) if fx == 0.0: return x if fb*fx < 0.0: a, fa = x, fx else: b, fb = x, fx return x # constante de los gases ideales R = 8.314472 # Constantes de van der Waals para el CO_2 a, b = 0.3640, 4.267e-5 # temperaturas (K) Tl = arange(100.0, 1000.1, 100.0) # presiones (atm) Pl = logspace(0.0, 2.0, 13) print("Ecuaci´ on de van der Waals comparada con la de los gases ideales") print("Gas: CO2, constantes de van der Waals (SI):") print("a = %e, b = %e\n" % (a, b)) for T in Tl: # barriendo el array de temperaturas print("\n %5s %8s %12s %12s %4s %12s %4s" % ("T ", "P ", "v ideal ", "v bisec ", "nfun", "v fsolve ", "nfun")) print(" (K) (atm) (L/mol) (L/mol) (L/mol)") for Patm in Pl: # barriendo el array de presiones P = 101325.0*Patm # volumen molar ideal videal = R*T/P # C´ alculo de v por el m´ etodo de bisecci´ on nfun = 0 vVWb = biseccion(fVW, 1.0e-6*videal, 10*videal) nfunb = nfun

Introducci´on a la F´ısica Computacional

16

Universidad de Oviedo

Departamento de F´ısica

# C´ alculo de v usando fsolve nfun = 0 vVWs = fsolve(fVW, videal) nfuns = nfun print(" %5d %8.2f %12.8f %12.8f %4d %12.8f %4d" % (T, Patm, 1.0e3*videal, 1.0e3*vVWb, nfunb, 1.0e3*vVWs, nfuns))

2.6. Ejemplos adicionales y de ´ındole pr´actica Otros m´etodos num´ericos para determinar ceros de funciones: m´etodo de Newton-Raphson, m´etodo del punto fijo, m´etodo de la secante, m´etodo de regula falsi. Funciones con comportamientos an´omalos. Extensi´on de estos m´etodos a casos multidimensionales.

Introducci´on a la F´ısica Computacional

17

Problema 3 Determinaci´on num´erica del per´ıodo de un p´endulo de amplitud cualquiera (Cuadratura num´erica en una dimensi´on) La obtenci´on de integrales de Riemann (o sea las denominadas ‘cuadraturas’ num´ericas—el nombre viene del problema de la cuadratura del c´ırculo, obtenci´on del a´ rea del c´ırculo) es un problema tambi´en b´asico en C´alculo Num´erico y F´ısica Computacional. Se necesita a menudo obtener una integral o un a´ rea entre curvas. Los m´etodos adecuados son realmente sumas ponderadas (s´olo que sumas de muchos valores de la funci´on). Los diversos algoritmos se diferencian t´ıpicamente en los factores de ponderaci´on de los valores de la funci´on en diversos puntos del intervalo de integraci´on. Tambi´en existen m´etodos que utilizan extrapolaciones al l´ımite. En este tema desarrollaremos un m´etodo sencillo para calcular integrales: un m´etodo con rect´angulos, y propondremos tambi´en otro m´etodo denominado de Simpson. Asimismo usaremos la funci´on correspondiente que aparece en scipy para cuadraturas unidimensionales.

3.1. Objetivos Estudiando el movimiento de un p´endulo de amplitud cualquiera se aborda el problema del c´alculo de la integral definida, para obtener su per´ıodo. Se plantea una soluci´on simple por el m´etodo de los rect´angulos y luego se aborda la soluci´on usando el m´etodo de Simpson y la funci´on quad de scipy.

3.2. Conocimientos previos Mec´anica cl´asica: oscilaciones. 18

Universidad de Oviedo

Departamento de F´ısica

Matem´aticas: c´alculo integral.

3.3. Planteamiento del problema El p´endulo es un s´olido r´ıgido que puede oscilar en un plano vertical alrededor de un eje de oscilaci´on perpendicular a ese plano. Cuando el s´olido r´ıgido se sustituye por una part´ıcula que oscila colgando de una cuerda el p´endulo se denomina ‘matem´atico’ o ‘simple’. El per´ıodo del p´endulo depende de la amplitud angular del movimiento de una manera compleja, la cual se va a analizar aqu´ı. Se plantea inicialmente la ecuaci´on de movimiento del p´endulo. Hag´amoslo para el p´endulo simple, aplicando la segunda ley de Newton a la fuerza tangencial:

111111111111111111 000000000000000000 000000000000000000 111111111111111111 000000000000000000 111111111111111111 000000000000000000 111111111111111111 000000000000000000 111111111111111111 000000000000000000 111111111111111111 000000000000000000 111111111111111111

θm

θ

F

t

Ft = −mg sen θ = mat = ml θ¨ El signo negativo es porque la fuerza tiende a disminuir el a´ ngulo. Operando, se obtiene:

θ¨ + ω02 sen θ = 0, siendo ω02 = g/l. Para el p´endulo f´ısico (se requiere estudiar din´amica del s´olido r´ıgido) se obtendr´ıa ω02 = (mgl)/I. En estas expresiones m es la masa del cuerpo que oscila, l la distancia del centro de masas del mismo al eje de oscilaci´on, I el momento de inercia del sistema respecto al eje de oscilaci´on y g la aceleraci´on de la gravedad del lugar.

3.4. Determinaci´on del per´ıodo de un p´endulo La expresi´on obtenida posee soluci´on anal´ıtica para amplitudes peque˜nas (y entonces θ ≈ sen θ ) pues se trata de un oscilador arm´onico. En ese caso se demuestra que el per´ıodo del p´endulo es T0 = 2π /ω0 . Sin embargo en el caso general el movimiento no es integrable. No obstante, se puede obtener una integral definida que nos d´e el per´ıodo del p´endulo. Puesto que en el movimiento participan s´olo fuerzas conservativas, la energ´ıa mec´anica se conserva: 0 + E p (θm ) = Ec (θ ) + E p (θ ), obteni´endose:

1 ω02 (1 − cos θm ) = θ˙ 2 + ω02 (1 − cos θ ), 2

Introducci´on a la F´ısica Computacional

19

Universidad de Oviedo

Departamento de F´ısica

donde θm es la amplitud del movimiento oscilatorio. Operando se llega f´acilmente a: √ p dθ θ˙ = = ω0 2 cos θ − cos θm . dt Por integraci´on de esta ecuaci´on entre los extremos (−θm , θm ) del movimiento se obtiene el per´ıodo: Z θm 1 1 √ T = T0 √ dθ , π 2 −θm cos θ − cos θm que es una integral impropia. Haciendo un cambio sencillo de variables se obtiene: 2 T = T0 π

Z

π 2

0

1 p dφ , 1 − k2 sen2 φ

(3.1)

con k = sen(θm /2). Esta integral se denomina ‘el´ıptica de primera especie completa’, la cual ya no es impropia excepto para k = 1 (´o sea para θ = π , que nos da infinito). Nuestro inter´es es realizar esta integral con valores diferentes de θm desde 0◦ hasta 180◦ , o π radianes, que es la amplitud m´axima posible, la cual nos da un per´ıodo infinito, ya que el p´endulo est´a en equilibrio inestable y no hay oscilaciones. Se trata de hacer un barrido de valores de amplitud, por ejemplo cada grado sexagesimal, para determinar el per´ıodo del p´endulo. La integral planteada se puede realizar num´ericamente de muy diversas maneras. Aqu´ı se exponen tres maneras diferentes.

3.5. M´etodos num´ericos de c´alculo de integrales definidas Queremos realizar el c´alculo num´erico de: S=

Z b

y f ( x)

f (x) dx.

a

Ese problema se denomina ‘integraci´on’ o ‘cuadratura num´erica’. El valor de S coincide con el a´ rea bajo la curva hasta el eje de abscisas entre a y b (v´ease la figura adjunta). Si dividimos el recinto de integraci´on en n partes (numeradas i = 0, 1, . . ., n − 1), cada una de anchura ∆x = h = (b − a)/n, podemos sustituir el a´ rea bajo la curva por el a´ rea del rect´angulo de esa anchura y de altura f (xi ), donde xi es el centro del intervalo i, o sea xi = a + (i + 12 )h. Esto da lugar al m´etodo de los rect´angulos centrados o m´etodo del medio punto: Z b a

n−1

f (x) dx ≈ h ∑ f (xi ).

S

a

b

x

y f ( x) h 0

a

1

i

xi

...

n −1

b

x

i=0

Introducci´on a la F´ısica Computacional

20

Universidad de Oviedo

Departamento de F´ısica

Otra metodolog´ıa t´ıpica se denomina m´etodo de Simpson, que consiste en dividir en intervalos de anchura h como antes pero tom´andolos de dos en dos (con tres abscisas, una a cada lado y otra en el justo medio) hacer pasar una par´abola por los tres puntos (interpolaci´on) y luego realizar la integral de la par´abola. Es necesario que n sea par. Sumando las a´ reas para cada dos intervalos se tiene la integral total, que resulta: Z b

f (x) dx ≈

a

siendo ahora xi = a + ih, y:

 1    4 ci =  2   1

h n ∑ ci f (xi), 3 i=0 i=0 i %2 = 1 i %2 = 0 i=n

Obviamente existen otros m´etodos m´as avanzados que est´an implementados en las bibliotecas de m´etodos num´ericos, tanto en scipy como en otros lenguajes. El m´odulo integrate incluye rutinas para realizar integrales definidas (cuadraturas). En 1D la funci´on m´as general dentro del m´odulo integrate es quad, que nos realiza la integraci´on de una funci´on de una variable entre dos l´ımites. Se invoca (hay m´as argumentos opcionales): quad(f, a, b, epsabs=valor) siendo f la funci´on a integrar, a y b los l´ımites de integraci´on (que pueden ser ±∞), y epsabs el error absoluto admisible en la integraci´on (por defecto 1.5e-8). Ejemplo 3.5.1. Cuadratura en una dimensi´on. La fuerza por unidad de altura en N m−1 sobre el m´astil de un barco de vela viene dada por:   z F(z) = 200 exp(−z/5), 2+z viniendo expresada la altura respecto al suelo en metros. Determine la fuerza total ejercida sobre el m´astil si tiene 10 m. Soluci´on from numpy import * # para exp() from scipy.integrate import quad def f (z): return 200.0*z/(2.0+z)*exp(-z/5.0) F = quad(f, 0.0, 10.0, epsabs=1.0e-6) print(F) # (resultado, error aproximado actual) (462.50466684953761, 4.2931842087556273e-09)

Introducci´on a la F´ısica Computacional

21

Universidad de Oviedo

Departamento de F´ısica

Ejercicio 3.5.1. La variaci´on de la aceleraci´on de la gravedad g con respecto a la altitud, y, viene dada por la expresi´on: R2 g= g0 (R + y)2 donde R = 6371 km es el radio de la tierra, y g0 = 9.81 m/s2 es la aceleraci´on de la gravedad al nivel del mar. La variaci´on de la energ´ıa potencial gravitatoria, ∆U , de un objeto que se eleva desde la tierra viene dada por: ∆U =

Z h

mgdy

0

Calcule la variaci´on de la energ´ıa potencial de un sat´elite de masa m = 500 kg que se eleva desde la superficie de la tierra hasta una altura de 800 km.

˜ de un programa para determinar el per´ıodo del 3.6. Diseno p´endulo Se trata fundamentalmente de preparar un programa principal que se encargue de la entrada de datos y de la salida de resultados, con un bucle que permite ir variando la amplitud del p´endulo, una funci´on que nos permite calcular el integrando de la integral que queremos calcular, y otra funci´on que tomando, como entrada la funci´on integrando, los l´ımites de integraci´on y el n´umero de intervalos usados, nos devuelva el valor de la integral. Esta estrategia es muy interesante ya que la rutina que nos calcula la integral es gen´erica y vale para calcular cualquier integral definida de una variable. Hemos de destacar la necesidad de hacer que h sea peque˜no, pues desde el punto de vista matem´atico la aproximaci´on se convierte en identidad cuando h → 0. El m´etodo de los rect´angulos es de primer orden (o sea el error disminuye seg´un h), mientras que el de Simpson es de tercer orden (disminuye seg´un h3 ). Sin embargo, desde el punto de vista del c´alculo num´erico, puesto que cuanto menor sea ∆x hay que realizar m´as sumas, y en cada suma se comete un error de redondeo, no es conveniente en la pr´actica disminuir mucho h. Se recomienda probar con valores crecientes del n´umero de intervalos usados, n = 10, 100, 1000, 10000, 100000, y compararlos. Los dos m´etodos citados deben compararse con quad de scipy. El c´odigo es el siguiente: # -*- coding: utf-8 -*# per´ ıodo de un p´ endulo de amplitud cualquiera from numpy import * from scipy.integrate import quad def f (x): # k es una variable global que tiene # que estar definida cuando se invoca f(x) Introducci´on a la F´ısica Computacional

22

Universidad de Oviedo

Departamento de F´ısica

return 1/sqrt(1.0-(k*sin(x))**2) def rectangulos (f, a, b, n=100): h = (b-a)/n x = linspace(a+h/2, b-h/2, n) return h*sum(f(x)) a, b = 0.0, pi/2.0 thm = arange(0.0, 180.0, 1.0) print("Salida generada por ’pendulo1.py’") print("Amplitud de un p´ endulo de amplitud no peque~ na") print("El m´ etodo de los rect´ angulos usa 100 intervalos") print("%8s %16s %16s" % ("th_max", "rect´ angulos", "quad")) for th in thm: k = sin(deg2rad(th/2)) # aqu´ ı se define k Ir = 2/pi*rectangulos(f, a, b) Iq = quad(f, a, b, epsabs=1.0e-15) print("%8.1f %16.12f %16.12f" % (th, Ir, 2/pi*Iq[0])) Tambi´en se deja como ejercicio para los estudiante representar gr´aficamente el per´ıodo frente a la amplitud con matplotlib.pyplot.

3.7. Ejemplos adicionales y de ´ındole pr´actica Otros m´etodos num´ericos para obtener integrales: m´etodo de Gauss y m´etodo de Montecarlo. Los errores de redondeo y de truncamiento.

Introducci´on a la F´ısica Computacional

23

Problema 4 Movimiento unidimensional: ca´ıda de un paracaidista (Integraci´on de ecuaciones diferenciales ordinarias) El movimiento, desde un punto de vista muy general, se define como la variaci´on de algo con el tiempo. En todos los campos de la ciencia se necesita determinar la evoluci´on de alguna magnitud con el tiempo (posici´on de una part´ıcula en mec´anica, la cantidad de un contaminante en un lugar, el valor de una acci´on en bolsa, etc.). El problema surge del desconocimiento de la definici´on de la funci´on que cumple la variable de inter´es con el tiempo, puesto que en general se conoce c´omo var´ıa con el tiempo a trav´es de ecuaciones diferenciales. El ejemplo m´as claro desde el punto de vista f´ısico son las ecuaciones de movimiento de una part´ıcula, que definen la aceleraci´on de la misma en funci´on de diversos par´ametros del problema usando las leyes de Newton. Se necesitan, pues, m´etodos precisos y estables para determinar unas variables de inter´es (variables dependientes) seg´un var´ıa el tiempo (variable dependiente). (En algunas ocasiones la variable dependiente no es el tiempo, sino una longitud u otra magnitud.) En este tema desarrollaremos el m´etodo de Euler (el m´as simple de todos) para resolver sistemas de ecuaciones diferenciales ordinarias (SEDO, aunque en ingl´es se usa la abreviatura ODE). Tambi´en enfocaremos el uso de otros m´etodos mejores y por supuesto el m´etodo que incorpora scipy. No abordaremos, por su dificultad y falta de tiempo, el estudio de los problemas diferenciales en derivadas parciales (en los que aparecen m´as de una variable independiente). El tema es bastante complejo y se aborda en cursos superiores. S´olo abordaremos problemas de valores iniciales (en los que se conocen las variables dependientes para un valor dado de la variable dependiente). No trataremos los denominados problemas de contorno.

24

Universidad de Oviedo

Departamento de F´ısica

4.1. Objetivos A partir del estudio del movimiento unidimensional de ca´ıda de un paracaidista se aborda el problema de la integraci´on de sistemas de ecuaciones diferenciales ordinarias (SEDO). Se plantea una soluci´on simple por el m´etodo de Euler (y su derivado de Euler-Cromer) y luego se aborda la soluci´on usando la funci´on odeint de scipy.

4.2. Conocimientos previos Mec´anica cl´asica: leyes de Newton. Matem´aticas: ecuaciones diferenciales ordinarias.

4.3. Planteamiento del problema Un paracaidista en su ca´ıda puede aproximarse por una part´ıcula que sufre dos fuerzas: la de la gravedad, aproximadamente constante, y otra de frenado por el aire, que depende de su velocidad. La fuerza de rozamiento puede ser aproximadamente proporcional al m´odulo de la velocidad (rozamiento de Stokes) o al cuadrado de ese m´odulo (rozamiento de Newton). Nuestro problema conlleva resolver la ecuaci´on diferencial de movimiento obtenida a partir de las leyes de Newton, obteniendo la posici´on y la velocidad para diversos valores del tiempo. Supongamos inicialmente que el avi´on est´a pr´acticamente est´atico. Si denominamos x a la distancia recorrida por el paracaidista desde que salt´o del avi´on (y creciente hacia abajo), suponiendo un rozamiento de Newton (el que sucede realmente en el r´egimen de velocidades de un paracaidista), por aplicaci´on de la 2a ley de Newton, se tiene: mg − β v2 = ma, donde m es la masa del paracaidista (con todo su equipo), β un coeficiente de rozamiento que depende de la aerodin´amica y de las propiedades del aire, v su velocidad (positiva hacia abajo) y a su aceleraci´on (tambi´en positiva hacia abajo). Denominando γ = β /m se tiene: a = g − γ v2 .

(4.1)

Esta es la ecuaci´on de movimiento, que es necesario integrar para determinar x(t), v(t) a partir de dos condiciones iniciales. Podemos tener condiciones iniciales x(t0 ) = x0 , v(t0 ) = v0 , o sea valores de posici´on y velocidad para un tiempo dado. En ese caso el problema a resolver se denomina ‘problema de valores iniciales’ (en ingl´es, IVP: initial value problem). Pero tambi´en podemos tener dos condiciones en diferentes instantes de tiempo. En ese caso el problema es algo m´as dif´ıcil de resolver (no lo abordaremos ahora) y se denomina problema de contorno (en ingl´es BVP: boundary value problem). Introducci´on a la F´ısica Computacional

25

Universidad de Oviedo

Departamento de F´ısica

Los m´etodos anal´ıticos s´olo son aplicables en ciertos casos (como el presente). En general no es posible obtener soluciones basadas en funciones conocidas y se deben utilizar m´etodos num´ericos (de los que existen muy variados).

4.4. El m´etodo de Euler Empecemos por el caso m´as simple, de integrar una ecuaci´on diferencial de primer orden, como: x′ =

dx = f (x,t), dt

(4.2)

donde t es la variable independiente y x la variable dependiente. Se parte de una condici´on inicial x(t = t0 ) = x0 . Una derivada temporal es un cociente incremental cuando el incremento temporal tiende a cero. Denominando x′ ≡ v por x simplificar: x1 − x0 , ∆t→0 ∆t

v0

v0 = l´ım

donde x0 es el valor de x para t = t0 y x1 es la misma variable para t = t1 = t0 + ∆t. Si no tomamos l´ımites y asumimos que ∆t es suficientemente peque˜no, entonces la igualdad se sustituye por una aproximaci´on. Despejando x1 se tiene:

(pendiente)

x1

real

x1

aproximada

x0 t0

t1

t

x1 ≈ x0 + v0 ∆t. Por tanto, conocida la derivada de la variable dependiente en un punto, se puede determinar aproximadamente la variable dependiente en un punto cercano. La ecuaci´on anterior expresa el m´etodo de Euler. En el m´etodo de Euler se demuestra que el error local es del orden de ∆t 2 , por lo que conviene usar incrementos temporales peque˜nos. Por otro lado, aunque ∆t sea peque˜no y la precisi´on se mantenga controlada, este m´etodo es un m´etodo no recomendable (tiene fundamentalmente un inter´es te´orico para problemas IVP), ya que es inestable (dependiendo de las condiciones del problema a resolver). Un m´etodo es estable si produce soluciones acotadas cuando la soluci´on exacta es acotada y es inestable cuando produce una soluci´on no acotada cuando la soluci´on exacta es acotada. Pasemos a preparar un algoritmo que exprese el m´etodo de Euler para integrar la ecuaci´on diferencial (4.2). Conocidas las condiciones iniciales, t0 y x0 , sustituyendo las aproximaciones

Introducci´on a la F´ısica Computacional

26

Universidad de Oviedo

Departamento de F´ısica

por igualdades, se hace la iteraci´on siguiente: xk+1 = xk + vk (xk ,tk ) ∆t, tk+1 = tk + ∆t. Aunque el error local es de segundo orden, el error global a lo largo de todo el c´alculo es de primer orden. Ejemplo 4.4.1. La velocidad de una part´ıcula en movimiento viene dada por la expresi´on: dx = −x3 − sint dt Use el m´etodo de Euler para resolver la ecuaci´on diferencial anterior con la condici´on inicial x(t = 0) = 0 hasta un tiempo final 10, usando 1000 pasos y dibuje el resultado. Soluci´on # -*- coding: utf-8 -*from numpy import * import matplotlib.pyplot as plt def f (x,t): return -x**3-sin(t) # par´ ametros del problema a, b, ni = 0.0, 10.0, 1000 tp = linspace(a, b, ni+1) # vector de tiempos # ------ AA -----h = tp[1]-tp[0] # intervalo de integraci´ on x = 0.0 xp = [] # lista con posiciones for t in tp: xp.append(x) x += h*f(x,t) # un paso del m´ etodo de Euler # ------ BB -----plt.plot(tp, xp) plt.xlabel("t") plt.ylabel("x(t)") plt.show() Ejercicio 4.4.1. Una barrera de protecci´on se sit´ua al final de un circuito con el objetivo de parar coches que han perdido el control. Esta barrera se ha dise˜nado de forma que la fuerza que la barrera aplica al coche viene dada en funci´on de la velocidad v y del desplazamiento x de la parte frontal del la barrera, seg´un la expresi´on: F = kv3 (x + 1)3 Introducci´on a la F´ısica Computacional

27

Universidad de Oviedo

Departamento de F´ısica

donde k = 30 kg s m−5 es una constante. Un coche con una masa m = 1500 kg, impacta contra la barrera de protecci´on a una velocidad de 90 km/h. Calcule y represente la velocidad del coche en funci´on de su posici´on para 0 m ≤ x ≤ 3 m, teniendo en cuenta que la variaci´on de la velocidad con el desplazamiento de la barrera viene dada por la expresi´on: dv −kv2 (x + 1)3 = dx m

4.5. Integraci´on de la ecuaci´on de movimiento en el caso unidimensional general Planteemos el uso del m´etodo de Euler para integrar las ecuaciones del movimiento, en un problema IVP. En este caso, puesto que la ecuaci´on diferencial del movimiento es de segundo orden, se tiene que convertir el problema en un sistema de ecuaciones diferenciales de primer orden. El problema: x′′ = a(x, x′ ,t) puede convertirse en un sistema, denominando v ≡ x′ , de la manera siguiente: v′ = a(x, v,t),

x′ = v.

De manera similar a la realizada con x, tambi´en: v1 − v0 , ∆t→0 ∆t

a0 = l´ım

donde v0 es el valor de la velocidad en t = t0 y v1 es el valor de la velocidad para t = t1 = t0 + ∆t. Si no tomamos l´ımites y asumimos que ∆t es suficientemente peque˜no, entonces la igualdad se sustituye por una aproximaci´on. Despejando v1 se tiene: v1 ≈ v0 + a0 ∆t, que se puede extender tambi´en a lo largo del tiempo. Podemos, por tanto, expresar el m´etodo de Euler para integrar el movimiento unidimensional de la siguiente manera. Conocidas las condiciones iniciales, t0, x0 y v0 , sustituyendo las aproximaciones por igualdades, se hace la iteraci´on siguiente: vk+1 = vk + a(xk , vk ,tk ) ∆t, xk+1 = xk + vk ∆t, tk+1 = tk + ∆t.

(4.3)

Existe una mejora (aunque peque˜na) del m´etodo citado, denominada m´etodo de Euler-Cromer, que consiste en utilizar en la determinaci´on de x no la velocidad al principio del intervalo, sino Introducci´on a la F´ısica Computacional

28

Universidad de Oviedo

Departamento de F´ısica

la que se acaba de determinar para el final del mismo. El m´etodo tiene mejores propiedades num´ericas que el de Euler (puesto que es semi-impl´ıcito), pero sigue siendo de primer orden: vk+1 = vk + a(xk , vk ,tk ) ∆t, xk+1 = xk + vk+1 ∆t, tk+1 = tk + ∆t.

4.6. M´etodos num´ericos para integrar sistemas de ecuaciones diferenciales ordinarias (SEDO) El m´etodo de Euler anteriormente descrito es el m´as simple de todos los m´etodos num´ericos dise˜nados para integrar num´ericamente ecuaciones diferenciales. Por ejemplo, su aplicaci´on a ecuaciones diferenciales en derivadas parciales da lugar a los m´etodos de diferencias finitas. No obstante existen m´etodos de orden superior que son m´as eficaces. Por ejemplo, el m´etodo del punto medio (o Runge-Kutta de orden dos, RK2), tiene dos pasos en cada ∆t. Aplicado a las ecuaciones del movimiento se tiene: 1

v∗ = vk + a(xk , vk ,tk ) ∆t, 2 1

x∗ = xk + vk ∆t, 2 1

t∗ = tk + ∆t, 2

vk+1 = vk + a(x∗ , v∗ ,t∗ ) ∆t, xk+1 = xk + v∗ ∆t, tk+1 = tk + ∆t. Globalmente este m´etodo es de orden dos (el de Euler es de orden uno), aunque en cada paso temporal necesita calcular dos veces la aceleraci´on. Otros m´etodos muy utilizados son los siguientes: Runge-Kutta, predictor-corrector, de extrapolaci´on, multipaso, y simpl´ecticos.

4.6.1. Integraci´on de SEDO en scipy La funci´on odeint del m´odulo scipy.integrate incluye algunos de estos m´etodos citados en el p´arrafo anterior, con control autom´atico del paso de integraci´on (para mantener un error estimado dado en la soluci´on). Se resuelve el problema de condiciones iniciales siguiente. Dado el sistema: dy0 = f0 (y0 , y1 , . . . ,t), dt dy1 = f1 (y0 , y1 , . . . ,t), dt ... dyn−1 = fn−1 (y0 , y1 , . . . ,t), dt Introducci´on a la F´ısica Computacional

29

Universidad de Oviedo

Departamento de F´ısica

y conocidas yi (t = 0) = y0i , determinar yi (t) para diversos valores de t. Se invoca de la siguiente manera: odeint(fun, y0, t, args) siendo fun(y, t) la funci´on que define el SEDO, y0 los valores iniciales y t un array con una secuencia de valores en los que se quiere la soluci´on. args son argumentos opcionales que se pasan (si se desea) a fun. La funci´on odeint tiene tambi´en otros argumentos opcionales. Por ejemplo rtol que es el error relativo buscado en el c´alculo, que por defecto es 1.5 × 10−8. Ejemplo 4.6.1. La velocidad de una part´ıcula en movimiento viene dada por la expresi´on: dx = −x3 − sint dt Use la rutina odeint para resolver la ecuaci´on diferencial anterior con la condici´on inicial x(t = 0) = 0 hasta un tiempo final 10, usando 1000 pasos y dibuje el resultado. Soluci´on El c´odigo fuente ser´ıa como el expuesto antes para el m´etodo de Euler pero incluyendo: from scipy.integrate import odeint en la parte superior y sustituyendo el fragmento de c´odigo entre las l´ıneas: # ------ AA -----# ------ BB -----por: xp = odeint(f, 0.0, tp) Ejemplo 4.6.2. Las ecuaciones diferenciales de Lorenz fueron obtenidas en el estudio de modelos climatol´ogicos. Son las siguientes: dx = σ (y − x), dt

dy = rx − y − xz, dt

dz = xy − bz, dt

donde σ , r y b son par´ametros constantes. Prepare un programa para integrar, usando 10000 puntos las ecuaciones anteriores con σ = 10, r = 28 y b = 8/3, en el rango t ∈ [0, 50] con las condiciones iniciales (x0 , y0 , z0 ) = (0, 1, 0). Dibuje el resultado. (¿Qu´e se puede deducir de la soluci´on dibujada?) Soluci´on # -*- coding: utf-8 -*from numpy import * import matplotlib.pyplot as plt from scipy.integrate import odeint Introducci´on a la F´ısica Computacional

30

Universidad de Oviedo

Departamento de F´ısica

def f (xyz, t): x, y, z = xyz[0], xyz[1], xyz[2] # para simplificar la notaci´ on return [sigma*(y-x), r*x-y-x*z, x*y-b*z] # par´ ametros del problema sigma, r, b = 10.0, 28.0, 8.0/3.0 ni = 10000 tp = linspace(0.0, 50.0, ni+1) # vector de tiempos xyz = odeint(f, [0.0, 1.0, 0.0], tp) plt.plot(tp, xyz[:,0], ’r’) plt.plot(tp, xyz[:,1], ’g’) plt.plot(tp, xyz[:,2], ’b’) plt.xlabel("t") plt.ylabel("x,y,z") plt.show() Salida gr´afica del ejemplo 50

50

x y

40

z

40

30

30

10

z

x,y,z

20

20

0

−10

10 −20

−30 0

10

20

30

40

t

50

0 −20

−15

−10

−5

0 x

5

10

15

20

En la figura de la izquierda se representan x, y, z frente al tiempo t, mientras que en la de la derecha se representa z frente a x.

˜ de un programa para integrar la ecuaci´on de 4.7. Diseno movimiento Una vez explicados los m´etodos de integraci´on de SEDO desde el punto de vista del c´alculo num´erico debemos plantear su codificaci´on en Python para el caso del problema del paracaidista, controlado por la ecuaci´on (4.1). Se trata fundamentalmente de preparar un programa principal que se encargue de la entrada de datos y de la salida de resultados, con un bucle que permita ir variando el tiempo desde un valor Introducci´on a la F´ısica Computacional

31

Universidad de Oviedo

Departamento de F´ısica

inicial t0 hasta un valor final t f , con incrementos Dt de salida de resultados. En cada iteraci´on se llama a la subrutina que nos realiza n pasos de integraci´on seg´un uno de los m´etodos citados, de tal manera que ∆t = Dt/n. Esto es debido a que ∆t debe ser, normalmente, m´as peque˜no que el intervalo de salida de resultados. Asimismo hemos de programar una funci´on que nos devuelva la aceleraci´on de la part´ıcula, teniendo presente que sus argumentos deben ser exclusivamente x, v y t (cualquier otro par´ametro debe ir como variable global). Esto es as´ı para no tener que modificar cada vez la subrutina de integraci´on del movimiento que, de esta manera, es gen´erica. El c´odigo en Python siguiente nos resuelve el problema mediante el m´etodo de Euler, seg´un las ecuaciones (4.3). Se deja como ejercicio resolver el problema mediante el m´etodo de Euler-Cromer y el RK2. # -*- coding: utf-8 -*# ca´ ıda de un paracaidista from numpy import * import matplotlib.pyplot as plt def a_paraca (x, v, t): return g-gamma*v**2 # g, gamma son variables globales def euler (acel, x, v, tv, n=10): nv, xv, vv = tv.size, zeros_like(tv), zeros_like(tv) xv[0], vv[0] = x, v for k in range(1, nv): t, Dt = tv[k-1], (tv[k]-tv[k-1])/float(n) for i in range(n): a = acel(x, v, t) t, v, x = t+Dt, v+a*Dt, x+v*Dt xv[k], vv[k] = x, v return (xv, vv) g = gamma =

9.80 # m/s^2 0.10 # m^-1

tv = linspace(0, 3, 30) x, v = 0, 0 # condiciones iniciales xv, vv = euler(a_paraca, x, v, tv) plt.xlabel("t (s)") plt.plot(tv, xv, tv, vv) plt.legend([’x (m)’, ’v (m/s)’], loc=’best’) plt.show()

Introducci´on a la F´ısica Computacional

32

Universidad de Oviedo

En la figura de al lado se tiene la soluci´on para los par´ametros expuestos en la misma. Utilizar odeint de la biblioteca scipy de Python conlleva preparar el c´alculo de la aceleraci´on y de la velocidad en la misma funci´on de c´alculo de las derivadas de las variables dependientes. Ser´ıa de la manera siguiente:

Departamento de F´ısica

paracaidista, =9 80m s2 , =0 10m
1 g

25

.

/

.

(m) (m s)

x v

/

20

15

10

5

0 0.0

0.5

1.0

1.5

(s)

2.0

2.5

3.0

t

def funder (y, t, g, gamma): ’’’ y es el vector [x, v] ’’’ x, v = y[0], y[1] return array([v, g-gamma*v**2]) La invocaci´on de odeint ser´ıa: yv1 = odeint(funder, array([x, v]), tv, args=(g,gamma)) xv1, vv1 = yv1[:,0], yv1[:,1] # para extraer posiciones y velocidades odeint es mejor que los otros m´etodos, ya que est´a compilada (es m´as r´apida). Utiliza m´etodos de orden cinco y superiores, y tambi´en un control autom´atico del paso para mantener un error (estimado) dado (por defecto 1.49012e-8 tanto el error absoluto como el relativo).

4.8. Ejemplos adicionales y de ´ındole pr´actica Otros m´etodos m´as precisos de integraci´on de ecuaciones diferenciales. Adaptaci´on del m´etodo de Euler para un control autom´atico del error.

Introducci´on a la F´ısica Computacional

33

Problema 5 C´alculo num´erico de centros de masa y de momentos de inercia (Cuadratura num´erica en dos dimensiones) En ocasiones se necesitan realizar cuadraturas en varias dimensiones (por ejemplo en el c´alculo de magnitudes extensivas). La existencia de varias dimensiones dificulta bastante el c´alculo. En primer lugar los l´ımites de integraci´on son curvas (en 2D), superficies (en 3D), hipersuperficies (en m´as dimensiones). Adem´as, puesto que se realizan sumas, el c´alculo en varias dimensiones multiplica de manera potencial el tiempo de c´alculo (y tambi´en la cantidad de memoria necesaria). En este tema desarrollaremos un m´etodo sencillo para calcular integrales en dos dimensiones: un m´etodo similar al de los rect´angulos en una dimensi´on.

5.1. Objetivos Desarrollar un programa que permita calcular el centro de masas y el momento de inercia de un cuerpo plano de forma cualquiera. En la resoluci´on del problema se aplican rutinas de mallado y algunos aspectos de vectorizaci´on (como la funci´on where de numpy).

5.2. Conocimientos previos Mec´anica cl´asica: centros de masa y momentos de inercia.

34

Universidad de Oviedo

Departamento de F´ısica

5.3. Planteamiento del problema El centro de masas y el momento de inercia son propiedades importantes de los cuerpos extensos, formados por un conjunto grande de part´ıculas, cuando se quiere analizar su comportamiento mec´anico como un todo. Por ello es interesante su c´alculo que, excepto en casos simples, suele ser bastante engorroso. Aqu´ı desarrollamos un programa de ordenador que permite hacer ese c´alculo para un cuerpo con forma plana. La generalizaci´on del programa al caso tridimensional es, por otro lado, sencilla. El caso bidimensional, que parece de limitada aplicaci´on, tiene en realidad mucha utilidad en resistencia de materiales, una parte de la Mec´anica de los s´olidos deformables. Un cuerpo formado por n part´ıculas individuales, cada una en una posici´on (xi , yi ) y con una masa mi , i = 1, 2, · · · , n, tiene como coordenadas de su centro de masas, G: xG =

1 m

n

∑ mixi ,

i=1

yG =

1 m

n

∑ miyi,

i=1

siendo m = ∑ni=1 mi la masa total, y su momento de inercia respecto al eje Oz es: n

Iz = ∑ mi (x2i + y2i ). i=1

Esas ecuaciones se pueden generalizar para cuerpos continuos, con densidad superficial conocida en cada punto σ (x, y), obteni´endose: xG yG

1 = σ (x, y) x dxdy, m Ω Z 1 = σ (x, y) y dxdy, Zm Ω Z

Iz =

Ω

σ (x, y) (x2 + y2 ) dxdy,

donde Ω indica el recinto de integraci´on, siendo la masa total: m=

Z

Ω

σ (x, y) dxdy.

Una vez conocidas estas magnitudes, por aplicaci´on del teorema de Steiner, es f´acil determinar el momento de inercia respecto a un eje paralelo al Oz que pase por el centro de masas, G: IG = Iz − m(x2G + y2G ).

5.4. Integraci´on num´erica en recintos multidimensionales Una vez obtenidas las expresiones matem´aticas que nos permiten conocer la masa, el centro de masas y el momento de inercia del cuerpo plano, se deben preparar para su uso en el ordenador. Debemos realizar las integrales indicadas anteriormente de un modo num´erico. Introducci´on a la F´ısica Computacional

35

Universidad de Oviedo

Departamento de F´ısica

Un m´etodo simple para realizar esa integraci´on es definir un rect´angulo de dimensiones tales que todo el recinto de integraci´on est´e incluido en aqu´el. Luego debemos dividir el rect´angulo en peque˜nos rectangulitos de anchura ∆x y altura ∆y, y suponer que en cada uno la densidad superficial es constante. y 11111 00000 00000 11111 00000 11111 C 00000 11111

yf A 11111 00000 00000 11111 00000 11111

k ∆y y0

f1 f2

∆x

i O

1111 0000 0000 1111 0000 1111 B 0000 1111

xf

x0

x

Si hay nx × ny rect´angulos, el centro de cada uno de ellos, (i, k), vendr´a definido por:     1 1 xi = x0 + i − ∆x, yk = y0 + k − ∆y, 2 2 donde x0 e y0 son los l´ımites izquierdo e inferior, respectivamente, del rect´angulo. El tama˜no de los rectangulitos vendr´a dado por: ∆x =

x f − x0 , nx

∆y =

y f − y0 , ny

siendo x f e y f los l´ımites derecho y superior, respectivamente, del rect´angulo. Ahora se puede sustituir de manera aproximada las integrales anteriores por sumas, teniendo en cuenta que si el centro del rect´angulo considerado est´a fuera del recinto debemos suponer σ = 0. Aqu´ı tenemos un peque˜no problema. El rect´angulo B de la figura est´a claramente dentro de la regi´on de inter´es, mientas que el C est´a claramente fuera. En el caso de un rect´angulo como el A que est´a en el borde de la regi´on de inter´es no tenemos m´as remedio que decidir si est´a fuera completamente o dentro completamente. Usaremos el criterio siguiente para ello: si el centro del rect´angulo est´a en la regi´on de inter´es entonces todo el rect´angulo estar´a dentro tambi´en. Esto significa que aproximaremos mediante poligonales con lados paralelos a Ox y a Oy las curvas l´ımite f1 (x) y f2 (x), como se ve en la figura siguiente. Esto implica un error, denominado de truncamiento, ya que estamos aproximando una funci´on continua por una aproximaci´on a trav´es de segmentos.

Introducci´on a la F´ısica Computacional

36

Universidad de Oviedo

Departamento de F´ısica y 1111 0000 0000 1111 0000 1111 0000 000001111 11111 00001111 00001111 1111 0000 00000 11111 0000 1111 0000 1111 0000 1111 00000 11111 00000 11111 0000 1111 0000 1111 0000 00000 11111 00000 11111 000001111 11111 00001111 00001111 0000 1111 00000 11111 00000 11111 00000 11111 0000 1111 0000 1111 0000 1111 00000 11111 00000 11111 00000 11111 0000 0000 00000 11111 0000 00000 11111 000001111 11111 00001111 1111 00001111 1111 0000 1111 00000 11111 00000 11111 00000 11111 0000 1111 0000 1111 0000 1111 00000 11111 00000 11111 00000 11111 0000 1111 0000 1111 0000 1111 00000 00000 11111 000001111 11111 00001111 0000 11111 00000 11111 00000 11111 0000 1111 0000 1111 00000 11111 000001111 11111 00001111 0000

f1 f2

O

x

Por tanto, las expresiones a codificar son: nx ny

m =

∑∑

σ (xi , yk ) ∆x∆y,

∑∑

σ (xi , yk ) xi ∆x∆y,

∑∑

σ (xi , yk ) yk ∆x∆y,

∑∑

σ (xi , yk )(x2i + y2k ) ∆x∆y,

i=1 k=1 nx ny

Sx =

i=1 k=1 nx ny

Sy =

i=1 k=1 nx ny

Iz =

i=1 k=1

y luego aplicaremos: xG =

Sx , m

yG =

Sy , m

IG = Iz − m(x2G + y2G ).

Del an´alisis de este m´etodo se deduce f´acilmente que es una generalizaci´on del m´etodo del medio punto al caso bidimensional. Una cuesti´on importante que se debe comentar es la dificultad de aumentar mucho el n´umero de rect´angulos (que es nx × xy ) si queremos que el error de redondeo (proporcional aproximadamente al n´umero de rectangulitos) no se dispare. Esto significa que al aumentar el n´umero de rectangulitos es posible que el error total vuelva a aumentar. La figura elegida como ejemplo es un cuarto de c´ırculo, y puesto que se conoce anal´ıticamente la soluci´on del problema planteado, es posible comparar la soluci´on anal´ıtica con la num´erica. Se deja como ejercicio el estudio de otras figuras geom´etricas con soluci´on anal´ıtica que se pueden encontrar en los libros de Mec´anica. El m´odulo integrate incluye una rutina para realizar integrales definidas en 2D que se denomina dblquad, la cual nos realiza la integraci´on de una funci´on de dos variables. Su uso es algo m´as complejo que quad pues, como ya hemos expuesto en el m´etodo simple antes dise˜nado, hacen falta dos tipos de l´ımites: para la abscisa, x, se deben dar dos l´ımites a y b, mientras que para la ordenada, y, hacen falta dos funciones, f1 (x) que nos da el l´ımite inferior, y f2 (x) que nos da el l´ımite superior. Introducci´on a la F´ısica Computacional

37

Universidad de Oviedo

Departamento de F´ısica

Se invoca (hay m´as argumentos opcionales): dblquad(f, a, b, f1, f2, epsabs=valor) siendo f la funci´on a integrar, a y b los l´ımites de integraci´on en abscisas (que pueden ser ±∞), f1 y f2 las funciones l´ımite en ordenadas, y epsabs el error absoluto admisible en la integraci´on (por defecto 1.5e-8). Las funciones f1 y f2 pueden definirse mediante def pero tambi´en se puede usar la instrucci´on lambda (como en el ejemplo expuesto m´as adelante). Ejemplo 5.4.1. Cuadratura en dos dimensiones. La temperatura de una placa rectangular se describe mediante la funci´on (las distancias en metros y la temperatura en kelvin): T (x, y) = 2xy + 2x − x2 − 2y2 + 40, La placa verifica x ∈ [0, 8] , y ∈ [0, 6]. Determine la temperatura media, sabiendo que: 1 T¯ = S

Z x2 Z y2 x1

T (x, y)dxdy.

y1

Soluci´on from numpy import * from scipy.integrate import dblquad def T (x, y): return 2*x*y + 2*x - x**2 - 2*y**2 + 40.0 S = 6.0*8.0 # superficie de la placa F = dblquad(T, 0.0, 8.0, lambda x: 0.0, lambda x: 6.0, epsabs=1.0e-6) print(F[0]/S) # resultado: 15.3333333333

˜ de un programa para el c´alculo de centros de 5.5. Diseno masa y momentos de inercia Se prepara una funci´on que tenga como argumentos los l´ımites inferior y superior, y el n´umero de puntos del mallado en cada direcci´on, y las funciones l´ımite superior e inferior. El programa principal se encarga de la entrada de datos y de la salida de resultados. El c´odigo, vectorizado en lo posible, es el siguiente:

Introducci´on a la F´ısica Computacional

38

Universidad de Oviedo

Departamento de F´ısica

# -*- coding: utf-8 -*# c´ alculo de centros de masa y momentos de inercia from numpy import * R = x0, y0, sig

1.0 # radio del c´ ırculo xf = 0.0, 1.0 # l´ ımites del rect´ angulo en Ox yf = 0.0, 1.0 # l´ ımites del rect´ angulo en Oy = 1.0 # densidad de masa

def f1 (x): return 0.0 def f2 (x): return sqrt(R**2-x**2) def cdmmominer (x0, xf, y0, yf, nx, ny, f1, f2): dx, dy = (xf-x0)/nx, (yf-y0)/ny x = linspace(x0+0.5*dx, xf-0.5*dx, nx) y = linspace(y0+0.5*dy, yf-0.5*dy, ny) yv, xv = meshgrid(y, x) s = where(yv