2.2: Mejorando el Programa

Volvamos al programa para calcular el potencial eléctrico, que discutimos en la parte anterior del tutorial, y mejorarlo aún más. Estas mejoras mostrarán algunas características más avanzadas de Python y Scipy que son buenas para conocer.

También haremos un cambio sustantivo en la física: en lugar de tratar las partículas como objetos puntuales, asumiremos que tienen un radio finito\(R\), con toda la carga concentrada en la superficie. De ahí que el potencial producido por una partícula de carga total\(q_{0}\) y posición\(x_{0}\) tendrá la forma

\[\phi(X) = \left\{\begin{array}{cl} \displaystyle \frac{q_0}{|X-x_0|},&\mathrm{if}\;\;|X-x_0| \ge R \\\displaystyle \frac{q_0}{R} &\mathrm{if}\;\;|X-x_0| < R. \end{array}\right.\]

Abra algunos archivos Python, y llámalo potentials2.py. Escribe en él lo siguiente:

from scipy import *
import matplotlib.pyplot as plt

## Return the potential at measurement points X, due to particles
## at positions xc and charges qc.  xc, qc, and X must be 1D arrays,
## with xc and qc of equal length.  The return value is an array
## of the same length as X, containing the potentials at each X point.
def potential(xc, qc, X, radius=5e-2):
    assert xc.ndim == qc.ndim == X.ndim == 1
    assert len(xc) == len(qc)
    assert radius > 0.

    phi = zeros(len(X))
    for j in range(len(xc)):
        dphi = qc[j] / abs(X - xc[j])
        dphi[abs(X - xc[j]) < radius] = qc[j] / radius
        phi += dphi
    return phi

## Plot the potential produced by N particles of charge 1, distributed
## randomly between x=-1 and x=1.
def potential_demo(N=20):
    X = linspace(-2.0, 2.0, 200)
    qc = ones(N)

    from scipy.stats import uniform 
    xc = uniform(loc=-1.0, scale=2.0).rvs(size=N)

    phi = potential(xc, qc, X)

    fig_label = 'Potential from ' + str(N) + ' particles'
    plt.plot(X, phi, 'ro', label=fig_label)
    plt.ylim(0, 1.25 * max(phi))
    plt.legend()
    plt.xlabel('r')
    plt.ylabel('phi')
    plt.show()

potential_demo(100)

Ahora vamos a pasar por los cambios que hemos hecho en el programa.

2.2.1 Parámetros de función opcionales

Como antes, definimos una función llamada potencial, cuyo trabajo es calcular el potencial eléctrico producido por una colección de partículas en 1D. Sin embargo, puede notar un cambio en el encabezado de la función:

def potential(xc, qc, X, radius=5e-2):

Hemos agregado un parámetro opcional, especificado como radius=5e-2. Un parámetro opcional es un parámetro que tiene un valor predeterminado. En este caso, el parámetro opcional se denomina radio, y su valor predeterminado es 5e-2 (lo que significa\(5\times10^{-2}\); también se puede escribir como 0.05, que es equivalente). Si llamas a la función omitiendo la última entrada, se asumirá que el valor es 0.05. Si proporciona un valor explícito para la última entrada, eso anula el valor predeterminado.

Si una llamada a función omite un parámetro no opcional (que como xc), eso es un error fatal: Python detendrá el programa con un mensaje de error.

2.2.2 Declaraciones de afirmación

En el cuerpo de la función, hemos agregado las siguientes tres líneas:

    assert xc.ndim == qc.ndim == X.ndim == 1
    assert len(xc) == len(qc)
    assert radius > 0.

La sentencia assert es una declaración especial de Python que verifica el valor de verdad de la siguiente expresión; si esa expresión es falsa, el programa se detendrá y se mostrará un mensaje de error informativo.

Aquí, utilizamos las declaraciones assert para verificar que

• xc, qc y X son todas matrices 1D (nota: el operador == Python comprueba la igualdad numérica)
• xc tiene la misma longitud que qc
• radio tiene un valor positivo (nota: 0. es Python short hand para el número 0.0).

Similar a escribir comentarios, agregar declaraciones de afirmación a su programa es una buena práctica de programación. Se utilizan para verificar que las suposiciones hechas por el resto del código (por ejemplo, que las matrices xc y qc tienen la misma longitud) se cumplen efectivamente. Esto asegura que si cometemos un error de programación (por ejemplo, suministrando matrices de tamaño incompatible como entradas), el problema surgirá lo antes posible, en lugar de dejar que el programa continúe ejecutándose y causando un error más sutil más adelante.

2.2.3 Rebanado avanzado

Dentro del bucle for, hemos cambiado la forma en que se calcula el potencial:

    for j in range(len(xc)):
        dphi = qc[j] / abs(X - xc[j])
        dphi[abs(X - xc[j]) < radius] = qc[j] / radius
        phi += dphi

Como se discutió anteriormente, ahora estamos considerando partículas de tamaño finito en lugar de partículas puntuales, por lo que el potencial es constante a distancias por debajo del radio de partícula. Esto se logra usando una técnica avanzada de corte en matriz.

Para cada partícula\(j\), el potencial se calcula en tres etapas:

• Calcule el potencial usando la fórmula regular\(q_j / |X-x_j|\), y guarde esos valores en una matriz, uno por cada valor de\(X\).
• Encuentre los índices de esa matriz que corresponden a los valores con\(|X - x_j| < R\), y sobrescriba esos elementos con el valor constante\(q_j/R\). Para encontrar los índices relevantes, hacemos uso de la siguiente característica de corte: si se proporciona una expresión de comparación como índice, se refiere a aquellos índices para los que la comparación es verdadera. En este caso, la expresión de comparación es abs (X-xC [j]) < radio, que se refiere a los índices de los\(X\) cuales están por debajo del radio mínimo. Estos índices son los que están en la matriz dphi que queremos sobrescribir.
• Agregar el resultado al potencial total.

Función Demo

Finalmente, tenemos una función de “demo” o (“demostración”) para hacer las parcelas adecuadas:

## Plot the potential produced by N particles of charge 1, distributed
## randomly between x=-1 and x=1.
def potential_demo(N=20):
    X = linspace(-2.0, 2.0, 200)
    qc = ones(N)

    from scipy.stats import uniform 
    xc = uniform(loc=-1.0, scale=2.0).rvs(size=N)

    phi = potential(xc, qc, X)

    fig_label = 'Potential from ' + str(N) + ' particles'
    plt.plot(X, phi, 'ro', label=fig_label)
    plt.ylim(0, 1.25 * max(phi))
    plt.legend()
    plt.xlabel('r')
    plt.ylabel('phi')
    plt.show()

potential_demo(100)

Mientras que nuestro programa anterior puso las cosas de plotting en “nivel superior”, aquí encapsulamos el código de plotting en una función potential_demo (). Esta función es llamada por la declaración de nivel superior potential_demo (100), que ocurre al final del programa.

Es útil hacer esto porque si, en el futuro, quieres que el programa demuestre algo más (por ejemplo, producir un tipo diferente de trama), no será necesario eliminar la función potential_demo (y arriesgarte a tener que reescribirla si cambias de opinión). En su lugar, puede escribir otra función de demostración y revisar esa única declaración de nivel superior para llamar a la nueva función de demostración en su lugar.

La función potential_demo proporciona otro ejemplo de uso de parámetros opcionales. Acepta un parámetro N=20, especificando el número de partículas a colocar. Cuando el programa se ejecuta, sin embargo, la función se invoca a través de la instrucción de nivel superior potential_demo (100), es decir, con una entrada real de la\(100\) cual anula el valor predeterminado de\(20\). Si la sentencia de nivel superior hubiera sido en su lugar potential_demo (), entonces se\(20\) usaría el valor predeterminado de.

Variables aleatorias de muestreo

La función demo genera N partículas con posiciones aleatorias. Esto se hace usando este código:

    from scipy.stats import uniform 
    xc = uniform(loc=-1.0, scale=2.0).rvs(size=N)

La primera línea importa una función llamada uniforme desde el módulo scipy.stats, que es un módulo que implementa distribuciones de números aleatorios. Como muestra este ejemplo, las declaraciones de importación no tienen que ser declaraciones de nivel superior. En algunos casos, podríamos optar por realizar una importación solo cuando se ejecuta una función en particular (generalmente, esto se hace si esa función es la única en el programa que depende de ese módulo).

La función uniforme devuelve un objeto que corresponde a una distribución uniforme particular. Uno de los métodos de este objeto, llamado rvs, genera una matriz de números aleatorios extraídos de esa distribución.

Trazado

Después de calcular el potencial total usando una llamada a la función potencial, lo trazamos:

   fig_label = 'Potential from ' + str(N) + ' particles'
   plt.plot(X, phi, 'ro', label=fig_label)

Para empezar, concéntrate en la segunda línea. Este es un uso un poco más sofisticado de la función plot de Matplotlib que el que tuvimos la última vez.

Los dos primeros argumentos, como antes, son las coordenadas x e y para la trama. El siguiente argumento, 'ro', especifica que queremos trazar usando círculos rojos, en lugar de usar líneas con el color predeterminado.

El cuarto argumento, label=fig_label, especifica algún texto con el que etiquetar la curva trazada. A menudo es útil asociar cada curva de una figura con una etiqueta (aunque, en este caso, la figura contiene solo una curva).

Esta forma de especificar una entrada de función, que tiene la forma FOO=BAR, es algo que no hemos visto previamente. Se basa en una característica conocida como argumentos de palabras clave. En este caso, label es la palabra clave (el nombre del parámetro que estamos especificando), y fig_label es el valor (que es un objeto string; lo discutiremos a continuación). Los argumentos de palabra clave permiten al llamador de una función especificar parámetros opcionales en cualquier orden. Por ejemplo,

   plt.plot(X, phi, 'ro', label=fig_label, linewidth=2)

es equivalente a

   plt.plot(X, phi, 'ro', linewidth=2, label=fig_label)

La lista completa de palabras clave para la función plot se da es su documentación.

Construyendo una cadena de etiqueta

A continuación, pasamos a la línea

   fig_label = 'Potential from ' + str(N) + ' particles'

que crea un objeto Python llamado fig_label, que se utiliza para etiquetar la curva. Este tipo de objeto se llama cadena de caracteres (o simplemente cadena para abreviar).

En el lado derecho de la declaración anterior, construimos el contenido de la cadena a partir de varias piezas. Esto se hace con el fin de obtener una cadena diferente para cada valor de N. El operador + “concatena” cadenas, uniendo las cadenas a su izquierda y derecha en una cadena más larga. En este caso, la cadena fig_label consta de las siguientes cadenas shofter, concatenadas juntas:

• Una cadena que contiene el texto 'Potencial de '.
• Una cadena que contiene el valor numérico de N, en forma de texto. Esto se calcula usando la función str, que convierte los números en sus representaciones de cadena correspondientes.
• Una cadena que contiene el texto 'partículas'.

El resto de la función potential_demo es relativamente autoexplicativa. La función ylim especifica los límites inferior y superior del\(y\) eje -axis de la gráfica (hay una función xlim similar, que no usamos). La declaración plt.legend () hace que la etiqueta de curva se muestre en una leyenda incluida en la trama. Finalmente, las funciones xlabel e ylabel agregan etiquetas de cadena a los\(y\) ejes\(x\) y.

Ejecutando el Programa

Ahora guarda tu programa potential2.py y ejecútalo. (Recordatorio: puedes hacerlo escribiendo python -i potential2.py desde la línea de comandos en GNU/Linux, o F5 en la GUI de Windows, o importar potencial2 desde la línea de comandos de Python). La cifra resultante se ve así:

Figura\(\PageIndex{1}\)