UNIVERSIDAD DE EL SALVADOR FACULTAD DE INGENIERÍA Y ARQUITECTURA ESCUELA DE INGENIERÍA DE SISTEMAS INFORMÁTICOS
TEMA: Integración en el Magnetismo. MATERIA: Análisis Numérico. CATEDRÁTICO: Lic. Guillermo Mejía Díaz. INTEGRANTES NOMBRE
CARNET
Morán Campos, Carlos Rolando
MC16006
Cerna Vásquez, Luis Roberto
CV16003
Rosales Alvarado, Joel Antonio
RA16004
Ciudad Universitaria, viernes 15 de junio de 2018.
TABLA DE CONTENIDO
INTRODUCCIÓN ....................................................................................................................... INTRODUCCIÓN ....................................................................................................................... 3 MARCO TEÓRICO DE REFERENCIA ...................................................................................... REFERENCIA ...................................................................................... 4 Reseña histórica del magnetismo....................................................................................... magnetismo. ...................................................................................... 4 Integración en magnetismo. ................................................................................................ magnetismo. ................................................................................................ 5 Ley de Biot-Savart. ............................................................................................................... Biot-Savart. ............................................................................................................... 5 Aplicaciones de la ley de Biot-Savart. ................................................................................ Biot-Savart. ................................................................................ 6 Integración numérica. .........................................................................................................17 numérica. .........................................................................................................17 Fórmulas de cuadratura de Newton-Cotes simples. simples. .........................................................18 Fórmulas de cuadratura de Newton-Cotes compuestas. .................................................19 compuestas. .................................................19 Integración numérica en Octave. .......................................................................................21 Octave. .......................................................................................21 DISEÑO DE LA SOLUCIÓN SOLUCIÓN .....................................................................................................24 MODELO COMPUTACIONAL COMPUTACIONAL . .................................................................................................25 FLUJOGRAMAS FLUJOGRAMAS ...................................................................................................................25 PSEUDO CÓDIGO CÓDIGO ................................................................................................................34 EVIDENCIA DE FUNCIONAMIENTO FUNCIONAMIENTO .......................................................................................41 .......................................................................................41 CONCLUSIONES CONCLUSIONES .....................................................................................................................43 RECOMENDACIONES RECOMENDACIONES .............................................................................................................44 BIBLIOGRAFÍA BIBLIOGRAFÍA ........................................................................................................................45 GLOSARIO GLOSARIO . ..............................................................................................................................46
2
INTRODUCCIÓN En análisis numérico, la integración numérica constituye una amplia gama de algoritmos para calcular el valor numérico de una integral definida y, por extensión, el término se usa a veces para describir algoritmos numéricos para resolver ecuaciones diferenciales. La integración numérica es una técnica que se puede usar para aproximar el valor de la integral de una función que no sea posible anti diferenciar (integrar). En el presente proyecto se desarrolla la solución de integración numérica aplicada al magnetismo, para ellos ocupamos la Ley de Biot Savart, “El campo magnético total generado por varias cargas en movimiento es la suma vectorial de los campos generados por las cargas individuales”, en palabras más simples, indica el campo
magnético creado por corrientes eléctricas estacionarias. Las aplicaciones de la Ley de Biot Savart, son campo magnético en las distintos tipos de forma, espira circular, corriente en forma general, solenoide, cascarón esférico con carga en rotación, cilindro hueco con carga en rotación. Para la solución de integrales se usa el método numérico de simpson compuesto para aproximación de integrales, consiste en aproximar los subintervalos de la función mediante polinomios de segundo grado. y con la cual resolvemos las aplicaciones de la Ley de Biot Savart. Para la codificación del método numérico para resolver integrales se usa Octave para su programación, y también para la creación de las interfaces gráficas usadas, se agrega que para la solución de aplicaciones que son integrales dobles, se usa el método dblquad, ya integrado en Octave.
3
MARCO TEÓRICO DE REFERENCIA
Reseña histórica del magnetismo. Los fenómenos magnéticos fueron observados por primera vez al menos hace 2500 años, con fragmentos de mineral de hierro magnetizado cerca de la antigua ciudad de Magnesia (hoy Manisa, en Turquía occidental). Esos trozos eran ejemplos de lo que ahora llamamos imanes permanentes; es probable que en la puerta del refrigerador de su hogar haya varios imanes permanentes. Vimos que los imanes permanentes ejercían fuerza uno sobre otro y sobre trozos de hierro que no estaban magnetizados. Se descubrió que cuando una varilla de hierro entraba en contacto con un imán natural, aquélla también se magnetiza, y si la varilla flotaba en agua o se suspendía de un hilo por su parte central, tendía a alinearse con la dirección norte-sur. La aguja de una brújula ordinaria no es más que un trozo de hierro magnetizado. El magnetismo, desde su origen y a lo largo de su historia logró mover los cimientos enteros de la humanidad. Este fenómeno encuentra su origen desde la Antigua Grecia, época en la cual, se observaron que ciertas piedras tenían la capacidad natural de atraer a objetos ferrosos. Tal como lo indican los hallazgos helénicos, el magnetismo consiste en el poder de atracción, o también, de repulsión, que tienen ciertos objetos sobre otros elementos caracterizados por la presencia de metales y componentes ferrosos, como por ejemplo, el níquel y el cobalto, al igual que la variedad de sus aleaciones. Mencionados elementos de energía magnética se le denominan como imánes, los cuales pueden ser naturales como los descritos durante la Antigua Grecia, así como ser artificiales, que logran la fuerza magnética tras atravesar diferentes procesos. Dentro de la historia griega, el filósofo Tales de Mileto se postula como el primer estudioso de este fenómeno de la naturaleza. Por su parte, en la literatura oriental, específicamente en el texto ‘Libro del amo del valle del diablo’, el cual data del siglo IV antes de Cristo, se mantiene una descripción de la magnetita por su po der para “atraer el hierro hacia sí o es atraída por este”.
Tras atravesar un largo periodo donde el magnetismo no se hallaba más allá del auxilio que ofrece a la navegación, llegando incluso a no abarcar otras áreas que estimularan un mayor uso de sus propiedades, el médico y físico, William Gilbert desarrolló en el año 1600, la investigación ‘De Magnete’, con la cual establece las bases primordiales
4
que catapultaron al magnetismo mediante la caracterización detallada de sus elementos, así como los tipos de imanes, y las posibilidades experimentales que tiene el fenómeno magnético. Gracias a Gilbert se conocen conceptos extraños con anterioridad, como lo fueron, los polos, los conductores, los aislantes, y la magnetización artificial. Tras los trabajos de Gilbert, otros científico como Hans Christian Orsted, André-Marie Ampère, Carl Friedrich Gauss, Michael Faraday, derivaron en el poder circulante del magnetismo en la naturaleza, así como su relación con el fenómeno de la luz, lo que antecedió a un cambio multidisciplinario de la manera en cómo el hombre lograba entender conceptos básicos como la energía, el espacio, y el tiempo, gracias a la teoría de la relatividad descrita por Albert Einstein.
Integración en magnetismo. El cálculo de los campos magnéticos de los alambres portadores de corriente se realiza usando la ley Biot-Savart. Esta ley se escribe necesariamente en términos de una integral, porque a diferencia de los campos eléctricos y gravitacionales, cuyas expresiones matemáticas implican cargas y masas puntuales locales, los campos magnéticos son causados por corrientes eléctricas extendidas y, como tales, requieren integrales.
Ley de Biot-Savart. Igual que en el campo eléctrico, hay un principio de superposición de campos magnéticos: “El campo magnético total generado por varias cargas en movimiento es la suma vectorial de los campos generados por las cargas individuales”
Comenzamos con el cálculo del campo magnético ocasionado por un segmento corto dl de un conductor que transporta corriente, como se ilustra en la siguiente figura:
El volumen del segmento es Adl, donde A es el área de la sección transversal del conductor, si hay n partículas con carga en movimiento por unidad de volumen, cada una con una carga q, la carga total dQ que se mueve en el segmento es:
5
=
La magnitud del campo eléctrico resultante es en cualquier punto a una distancia r es:
En forma vectorial:
Las dos ecuaciones anteriores son las que constituyen la Ley de Biot y Savart, pero eso es para un diferencial de campo magnético, para encontrar el campo magnético B total se integra la ecuación, simbólicamente es la siguiente:
Aplicaciones de la ley de Biot-Savart. Espira Circular. Si se mira en el interior de un timbre para puerta, un transformador, un motor eléctrico o un electroimán, se encontrarán bobinas de alambre con gran número de vueltas, espaciadas tan estrechamente que cada vuelta está muy cerca de formar una espira plana circular. En tales bobinas se utiliza una corriente para establecer un campo magnético. Por ello, es conveniente obtener una expresión para el campo magnético que produce una sola espira conductora circular portadora de corriente, o para las N espiras circulares estrechamente espaciadas que forman la bobina. Comencemos con el ejemplo simple de un bucle circular que lleva corriente de radio a ubicado en el plano xy. Con el centro del círculo en el origen de un sistema de coordenadas cilíndricas, los tres componentes cilíndricos del campo magnético son:
() = ∫ ( + −2() + )3⁄2 = 0 (() − ) = − ∫ ( + −2() + )3⁄2 Tecleamos los componentes que nos aparecen así:
6
Donde hemos tomado el radio del ciclo para que sea la unidad e ignoramos las constantes delante de la integral. Podríamos haber usado NIntegrate en su lugar; sin embargo, como resultado, Mathematica puede hacer las integrales analíticamente, en términos de las funciones elípticas. Es instructivo representar un gráfico de las líneas de campo del campo magnético. Para hacerlo, primero tenemos que cargar el paquete de gráficos especial que dibuja las líneas de campo. Esto se hace escribiendo primero << Graphics'PlotField ' Y luego el comando:
Corriente en forma general. El ejemplo de un bucle circular fue muy simple. Teníamos las integrales a nuestra disposición, y todo lo que teníamos que hacer era evaluarlas. Mathematica puede hacer mucho mejor que eso. Pero requiere que derivemos una fórmula general para el campo magnético de un filamento portador de corriente cuya forma es bastante arbitraria. Usamos coordenadas cartesianas para simplificar. Supongamos que el filamento se describe en coordenadas cartesianas mediante una ecuación paramétrica de la forma El resultado es el diagrama que se muestra en la Figura --. La figura muestra solo una sección transversal de las líneas de campo cortadas por un plano perpendicular a la espira y pasando por su centro. La sección transversal de la espira se muestra como dos espacios vacíos a mitad de camino a la derecha e izquierda de la vertical del medio.
7
Figura n. Las líneas del campo magnético de un bucle circular. El diagrama muestra solo las líneas tal como aparecen en un plano perpendicular al bucle y que pasan por su centro.
Solenoide. Con las fórmulas generales para los tres componentes de B a nuestra disposición, podemos calcular los campos de corrientes de forma específica. Una caso de interés es la hélice, correspondiente a un solenoide. Vamos a escribir:
Esto describe un solenoide de radio 1, el principio y el final de cuyo filamento ocurren en (cos (− 5), sin (− 5), 0.05 × (− 5)) [o (0,284, − 0,96, − 0,25)] y (0,284, − 0,96, 0,25),
respectivamente. Así pues, la longitud del solenoide es 0,5, y el espaciamiento entre los bobinados consecutivos es 0.05 × π = 0,314. Por lo tanto, hay por lo menos de dos
bobinas en el solenoide. Tracemos los componentes x y z del campo magnético fuera del solenoide. Elegimos un plano paralelo al plano XY que cruza el eje z en z = 1, y la ploteamos Bx y Bz en función de x, es decir, para los puntos de campo a lo largo del eje x. Dado que tenemos una expresión general para el campo, y las funciones que describen el solenoide están todas definidas anteriormente, todo lo que tenemos que hacer es trazar las funciones.
8
Figura 3.12. El componente z del campo magnético dentro de un solenoide de radio 1. El diagrama superior izquierdo es el campo del solenoide cuando su longitud L es 0,5. Moviendo hacia la derecha, vemos Bz para L = 1, L = 2, L = 5, L = 8, y L = 10, respectivamente.
Produce el ploteo a la izquierda de la figura 3,11 y
Produce el ploteo a la derecha. Estos ploteos demuestran el comportamiento de los componentes esperados sobre fundamentos físicos. Por ejemplo, se espera que Bx sea máximo cuando el punto está justo encima del alambre (casi una sola espira) y Bz para ser máximo en el centro. Ahora concentrémonos en los puntos cercanos al plano XY. De hecho, vamos a encontrar los componentes del campo en ese plano. Esto nos ayudará a investigar el límite del solenoide largo del campo, que se calcula en los cursos de física elemental usando la ley de circuitos de Ampere.El lector debe recordar que para tal solenoide, el campo interior es constante y enteramente en la dirección z, y fuera de él es cero. Es la constancia del campo que nos gustaría investigar a medida que aumenta la longitud del solenoide. Empezamos con el solenoide corto de longitud L = 0,5 como arriba. Una tentativa en calcular Bx o By producirá una queja de Mathematica sobre la integral siendo oscilatoria y convergiendo demasiado lentamente, pero la respuesta dada hacia fuera es casi cero.
9
Este campo magnético "constante" puede ser sorprendente, ya que la longitud del solenoide es una única fracción de su radio, y difícilmente podría ser llamado "largo". Sin embargo, un solenoide muy corto puede ser aproximado por una espira circular; y para una espira circular en el plano XY, el campo es completamente en la dirección z [ver ecuación (3,24)]. Sin embargo, el componente z del campo magnético puede calcularse numéricamente y trazarse. El comando
produce el diagrama que se muestra en la esquina superior izquierda de la figura 3.12. Esta figura muestra claramente que Bz tiene alguna variación notable a medida que se mueve lejos del eje del solenoide hacia su superficie lateral. El resto de la figura 3.12 muestra lo que sucede cuando la longitud del solenoide aumenta. Nuestra expectativa — basada en nuestra experiencia en cursos introductorios de física — es que el campo debe mostrar menos y menos sensibilidad a x a medida que la longitud del solenoide aumenta.Y esta expectativa se lleva a cabo en la figura. Salvo la fluctuación en el último diagrama — causada por fallas en los cálculos numéricos — incluso para longitudes moderadas de 8 y 10 (sólo 4 y 5 veces el diámetro del solenoide), Bz parece ser bastante constante para una buena fracción de la distancia desde el eje a la superficie lateral.El tratamiento elemental de un solenoide largo nos dice que el componente z (el único componente nonevanescente) del campo magnético debe seguir siendo constante a lo largo de la dirección z también. Así, Si trazamos Bz en función de z, debe permanecer bastante constante dentro del solenoide. La insensibilidad del campo a ambos x y z se puede exhibir mejor en un diagrama tridimensional de Bz en función de x y z donde el campo debe aparecer como una hoja plana (a excepción de los valores de x y z que se acercan a los bordes del solenoide). Un ploteo tridimensional puede ser producido escribiendo
El resultado es la figura 3.13. Observe que, aparte de la variación algo fuerte en la superficie lateral (que corresponde a x cerca de 1), Bz es constante relativamente en ambas direcciones, de x y z.
10
Cascarón esférico con carga en rotación. Cualquier carga en movimiento produce un campo magnético. Este movimiento puede ser en forma de una corriente generada por una batería o causada por un agente mecánico que actúa sobre cargas por lo demás estáticas. La fórmula general para este último tipo de campo magnético es:
donde dq (r') es el elemento de carga en r' y v (r ') su velocidad. Estamos interesados en el campo magnético generado por una capa esférica de radio uniformemente cargada que gira alrededor de uno de sus diámetros con una velocidad angul ar de ω. Debido a la geometría esférica de la fuente, usamos coordenadas esféricas para las variables de integración. Entonces nosotros tenemos,
Para las manipulaciones de Mathematica, es conveniente expresar todos los vectores unitarios en términos de vectores unitarios cartesianos con coeficientes esféricos,
FIGURA 3.13. La componente z del campo magnético dentro de un solenoide de radio 1, dibujado en función de x (distancia desde el eje) y z (distancia desde el centro en el eje). coordenadas Estos son:
11
Para el punto de campo y expresiones similares (con primo en las coordenadas) para el punto fuente. La razón para hacer esto es que Mathematica trata los vectores como teniendo componentes cartesianos. Por ejemplo, agrega componentes para encontrar los componentes de la suma, algo que generalmente no está permitido en las llamadas "coordenadas curvilíneas". Debido a la simetría del problema, no esperamos que el campo magnético dependa de φ de el punto de campo. Entonces, para simplificar, tomamos φ para que sea cero, es decir, posicionamos el punto de campo en el plano
xz. Ahora podemos escribir el código relevante para Mathematica. Primero, definimos los vectores de unidades esféricas en el punto de campo (con φ = 0):
donde usamos t para θ yp para φ. 7 Para el punto de origen usamos 1 en lugar de
primo:
También necesitamos los vectores unitarios cartesianos para encontrar los componentes cartesianos del campo:
La velocidad y todo el integrando ahora se pueden escribir en:
Para diferentes componentes correspondientes del integrando:
del
campo,
12
necesitamos
los
componentes
Finalmente, integramos numéricamente estos integrandos para obtener los componentes del campo:
Donde establecemos todas las constantes no geométricas multiplicando la integral igual a la unidad. Si escribimos Bx [0, 0, 1] — para el componente x del campo magnético en el centro de una esfera de radio 1 — conseguimos 9,67208 × 10 ^-18, indicando que Bx = 0 en el centro. Del mismo modo, conseguimos cero para By. Si cambiamos el punto de campo, seguimos obteniendo números pequeños para Bx (y By) siempre y cuando nos limitemos a los puntos interiores de la esfera. Por ejemplo, Bx [0.5, 1.5, 1] produce 4,80235 × 10^-8 y Bx [0.25, 1, 1] produce 2,3899 × 10^-10. Los componentes x e y del campo magnético de la cáscara esférica giratoria se desvanecen en el interior. El componente z es distinto de cero, sin embargo, y un ploteo de Bz en función de r para θ = 0 se muestra en la figura 3,14 donde r puede variar entre 0 y 2. El gráfico de la
izquierda, se ha producido escribiendo
muestra el resultado del cálculo en detalle indicando que el valor calculado de BZ está muy cerca de 8,37758, pero no exactamente el mismo para todos los puntos dentro. El diagrama de la derecha muestra todos los valores de BZ, incluidos los de afuera. se obtiene escribiendo
13
Figura 3,14. El componente z del campo magnético dentro de una cáscara esférica uniformemente cargada y uniformemente girando de radio 1, dibujada en función de r para θ = 0.
La línea horizontal para r entre 0 y 1 indica que Bz es constante en el interior — al menos para los puntos del eje polar (θ = 0). Una muestra de valores de
BZ para otros puntos dentro indica que es constante. Por ejemplo,
Todos producen 8.33758. Para los puntos fuera, el campo cae a cero a medida que el punto de campo se mueve más y más lejos de la esfera. Este comportamiento es evidente en el diagrama de la derecha de la figura 3.14.
Cilindro hueco con carga en rotación. Como último ejemplo de magnetismo, consideramos una concha cilíndrica hueca de radio uniformemente y longitud L que gira con una velocidad angular constante de ω.
Para el sistema de coordenadas cilíndrico usamos la siguiente ecuación,
y expresamos todos los vectores unitarios en términos de vectores unitarios cartesianos con coeficientes en coordenadas cilíndricas. Las coordenadas cartesianas y cilíndricas comparten el mismo êz. Los otros dos vectores unitarios se relacionan de la siguiente manera:
14
Como en el caso de la esfera giratoria, la simetría azimutal de la carcasa cilíndrica evita que el campo magnético dependa de φ del punto de campo. Entonces, una vez más, tomamos φ para que sea cero. El código relevante para Mathematica comienza con la definición de los vectores de unidades cilíndricas en el punto de campo (con φ = 0):
Nuevamente para el punto de origen usamos 1 en lugar de primo:
y define los vectores unitarios cartesianos como de costumbre:
La velocidad y todo el integrando ahora se pueden escribir en:
Para diferentes componentes correspondientes del integrando:
del
campo,
necesitamos
los
componentes
Finalmente, como antes, establecemos todas las constantes no geométricas que multiplican la integral igual a la unidad, e integramos numéricamente los integrandos para obtener los componentes del campo:
15
Estamos interesados en el comportamiento de B z a medida que el punto de campo se mueve en el plano xy. Esto podría hacerse usando un comando como
para un cilindro de unidad de radio y longitud L = 10. Sin embargo, debido a los pequeños valores de la integral exterior, Mathematica requeriría mucho tiempo para hacer la trama. La alternativa juiciosa es hacer una tabla con los valores seleccionados de r, y, si se desea, trazar la tabla usando ListPlot.
TABLA 3.1. Valores (c[i]) del campo magnético de un cilindro corto y un cilindro largo (b[i]). El primer cilindro tiene radio de unidad, longitud 0.5, y su campo magnético se denota por c [i]. Para obtener el quinto valor, por ejemplo, de Bz correspondiente a r [5] (que es 0.9) -para este cilindro, uno escribe: c[5]=Bz[0.9,0,0.5,1] El segundo cilindro tiene un radio de unidad, longitud 50, y su campo magnético se denota por b [i]. Para obtener el octavo valor de B z -correspondiente a r [8] (que es 1.01): para este cilindro, uno escribe: b[8]=Bz[1.01,0,50,1] En entradas separadas, asignamos valores a r. Por ejemplo, r[1]=0; r[2]=0.2; r[3]=0.5; r[4]=0.8; y entradas similares para r [5] a r [14]. Para construir la tabla de valores, tecleamos
16
Table[{r[i], c[i], b[i]}, {i, 1, 15}] // MatrixForm y Mathematica produce una tabla de tres columnas, que hemos formato de seis columnas) en la Tabla 3.1. Está claro que este características de un solenoide largo, es decir, campo magnético interior (la mitad izquierda de la mesa) y casi cero campo magnético mitad derecha de la mesa).
reproducido (en último tiene las constante en el en el exterior (la
Integración numérica. Introducción. Fórmula de cuadratura. Las fórmulas de integración numérica o de cuadratura son de la forma:
donde x0 x 1 0
1
x N (nodos) son N +1 puntos distintos pertenecientes al intervalo [a b] y N (pesos) son números reales.
Fórmula interpolatoria. Si P N es el polinomio que interpola a f en los puntos distintos x0 x 1 y
decimos que la fórmula de cuadratura es de tipo interpolatorio.
17
x N a b
Grado de precisión. Una fórmula de cuadratura tiene grado de precisión r si es exacta para
pero no es exacta para :
Fórmulas de cuadratura de Newton-Cotes simples. Son fórmulas de cuadratura de tipo interpolatorio, eligiendo los puntos de interpolación (nodos de la fórmula) igualmente separados de una de las dos formas siguientes: Fórmulas cerradas: los límites de integración a y b son nodos de la fórmula. ●
●
Fórmulas abiertas: ninguno de los límites de integración es nodo de la fórmula.
Fórmula de Simpson. La fórmula de Simpson simple es:
Usar la fórmula de Simpson para integrar una función en un intervalo, equivale a sustituir, dentro de la integral, la función a integrar por el polinomio de interpolación de
18
grado dos, que pasa por los puntos de la función de los extremos y el punto medio del intervalo. Es decir, sustituimos la función por una parábola e integramos.
Fórmulas de cuadratura de Newton-Cotes compuestas. Estas se obtienen dividiendo el intervalo ab en n subintervalos y aplicando a cada uno de estos subintervalos una fórmula de cuadratura sencilla. Sea:
19
Fórmula de Simpson compuesta:
20
Integración numérica en Octave. Octave viene con varias funciones incorporadas para calcular la integral de una función numéricamente (denominadas cuadratura). Todas estas funciones resuelven problemas de integración de 1 dimensión.
Funciones de una Variable. Octave soporta cinco algoritmos diferentes para calcular la integral de una función f sobre el intervalo de ‘a’ a ‘b’. Estos son: quad
Integración numérica basada en la cuadratura de Gauss. quadv
Integración numérica usando una regla de Simpson vectorizada adaptativa. quadl
Integración numérica mediante una regla de Lobatto adaptativa. quadgk
Integración numérica mediante una regla de Gauss-Konrod adaptativa. quadcc
Integración numérica usando reglas adaptativas de Clenshaw-Curtis. trapz, cumtrapz
Integración numérica de datos mediante el método trapezoidal. El mejor algoritmo de cuadratura para usar depende del integrando. Si tiene datos empíricos, en lugar de una función, la elección es trapz o cumtrapz. Si no está seguro de las características del integrando, quadcc será el más robusto, ya que puede manejar discontinuidades, singularidades, funciones oscilantes e intervalos infinitos. Cuando el integrando es uniforme quadgk puede ser el más rápido de los algoritmos. Función
Característica
quad
Baja exactitud con integrandos no uniformes.
21
quadv
Mediana exactitud con integrandos uniformes.
quadl
Mediana exactitud con integrandos uniformes. Más lento que quadgk.
quadgk
Mediana exactitud (1e-6 – 1e-9) con integrandos uniformes. Maneja las funciones oscilantes y los límites infinitos.
quadcc
Baja a alta precisión con integrands no uniformes/uniformes. Maneja funciones oscilantes, singularidades y límites infinitos.
Aquí hay un ejemplo del uso de la función quad para integrar la función: f(x) = x * sin (1/x) * sqrt (abs (1 - x))
De x=0 a x=3. El primer paso es definir la función: function y = x endfunction
.*
y sin
= .*
(1./x)
sqrt
f (abs
(1
-
(x) x));
El segundo paso es llamar a quad con los límites de la integración: [q,
ier,
nfun,
err]
=
quad
("f",
0,
3) 1.9819 1 5061
1.1522e-07
Funciones de múltiples variables. Octave no tiene funciones integradas para calcular la integral de funciones de múltiples variables directamente. Sin embargo, es posible calcular la integral de una función de variables múltiples usando las funciones existentes para integrales unidimensionales. Para ilustrar cómo se puede realizar la integración, integraremos la función f (x, y) = sin (pi * x * y) * sqrt (x * y)
para x y y entre 0 y 1.
22
El primer enfoque crea una función que integra f con respecto a x , y luego integra esa función con respecto a y . Debido a que quad está escrito en Fortran, no se puede llamar recursivamente. Esto significa que quad no puede integrar una función que llama quad y, por lo tanto, no puede usarse para realizar la integración doble. Sin embargo, se puede usar cualquiera de los otros integradores, que es lo que demuestra el siguiente código. función q = g (y) q = unos (tamaño (y)); para i = 1: longitud (y) f = @ (x) sin (pi * x. * y (i)). * sqrt (x. * y (i)); q (i) = quadgk (f, 0, 1); endfor función final I = quadgk ("g", 0, 1) 0.30022
El proceso anterior se puede simplificar con las funciones dblquad y triplequad para integrales en dos y tres variables. Por ejemplo: I = dblquad (@ (x, y) sin (pi * x. * Y). * Sqrt (x. * Y), 0, 1, 0, 1) 0.30022
: dblquad ( f , xa , xb , ya , yb ) : dblquad ( f , xa , xb , ya , yb , tol ) : dblquad ( f , xa , xb , ya , yb , tol , quadf ) : dblquad ( f , xa , xb , ya , yb , tol , quadf , ...) Evalúa numéricamente la integral doble de f . f es un identificador de función, función en línea o cadena que contiene el nombre de la función a evaluar. La función f debe tener la forma z = f (x, y) donde x es un vector e y es un escalar. Debería devolver un vector de la misma longitud y orientación que x . xa , ya y xb , yb son los límites inferior y superior de integración para xey respectivamente. El integrador subyacente determina si se aceptan límites infinitos. El argumento opcional tol define la tolerancia absoluta utilizada para integrar cada subintegral. El valor predeterminado es 1e-6. El argumento opcional quadf especifica qué función integradora subyacente usar. Cualquier opción, pero quad está disponible y el valor predeterminado es quadcc . Argumentos adicionales, se pasan directamente a f . Para usar el valor predeterminado para tol o quadf, uno puede pasar ':' o una matriz vacía ([]).
23
DISEÑO DE LA SOLUCIÓN Se realizará un software en octave para hacer los cálculos de las integrales que surgen de distintas aplicaciones de la Ley de Biot-Savart . Entre dichas aplicaciones se resolverán los siguientes casos: ● Corriente en una espira circular ● Corriente con forma general ● Solenoide ● Carcasa esférica cargada rotatoria ● Cilindro hueco cargado rotatorio Se creará una interfaz gráfica donde el usuario introduzca los parámetros requeridos para calcular el campo magnético. Por ejemplo en el caso de la espira circular hay que conocer el radio de la espira, la distancia r entre la espira y el punto donde se calculará el campo magnético, y la corriente eléctrica; cuando ya se han introducido los datos se ocupará la integral de la Ley de Biot-Savarta (para el caso de un cable conductor) para hacer los cálculos, y es la siguiente:
Para la solución numérica de la integral se usará el método de Simpsom, trapecio, o Simpsom compuesto, dejando la elección del método al usuario y así puede comparar cual de los métodos es el mejor. Además presentaremos resultados en forma gráfica:
Además se presentará una tabla con ciertos puntos, pero el usuario podrá interpolar algún valor que quiera calcular, esto se hará a través de trazadores cúbicos.
24
MODELO COMPUTACIONAL FLUJOGRAMAS Circular Loop
25
Current with general shape
26
27
Solenoide
28
29
Cylinder
30
31
Spherical
32
33
PSEUDO CÓDIGO Circular Loop { rho, z, limi, limf, N son variables con punto flotante } 1. Leer rho, z 2. Hacer Frho ← cos(t)/(((rho)^2+1 -(2*rho*cos(t))+(z)^2)^(3/2)) 3. Hacer Fz ← ((rho*cos(t))-1)/(((rho)^2+1-(2*rho*cos(t))+(z)^2)^(3/2)) 4. Hacer H ← (B -A)/N 5. Hacer XI0 ← Frho(A) + Frho(B) 6. Hacer XZ0 ← Fz( A) + Fz(B) 7. Hacer XI1 ← 0.0 8. Hacer XZ1 ← 0.0 9. Hacer XI2 ← 0.0 10. Hacer XZ2 ←0.0 11. Hacer NN ← N – 1 12. Hacer I ← 1
13. Repetir con I desde 1 hasta NN Hacer X ← A +( I * H) 14. 15. Si I MOD 2 es igual a 0 Hacer XI2 ← XI2 + Frho(X) 16. 17. Hacer XZ2←XZ2+Fz(X) 18. Sino Hacer XI1 ← XI1 + Frho(X) 19. Hacer XZ1 ← XZ1 + Fz(X) 20. 21. Fin Si Hacer I ← I +1 22. 23. Fin de ciclo paso 13 24. Hacer XI ← (XI0 + 2.0 * XI2 + 4.0 * XI1) * H / 3.0 25. Hacer XZ ← (XZ0 + 2.0 * XZ2 + 4.0 * XZ1) * H / 3.0 26. Hacer Bp ← XI 27. Hacer Bz ← XZ
34
Current with general shape { cx, cy, cz, I, limi, lims, ene son variables con punto flotante } { f, g, h son funciones } 1. Hacer ex ← [1, 0, 0] 2. Hacer ey ← [0, 1, 0] 3. Hacer ez ← [0, 0, 1] 4. Hacer r ← [x, y, z] 5. Hacer rp ← [f(t), g(t), h(t)] 6. Hacer x ← cx 7. Hacer y ← cy 8. Hacer z ← cz 9. Hacer Km ← 1/(4*pi) 10. Hacer intB_N ← cross( diff( rp(t), t), r(x, y, z) - rp(t) ) 11. Hacer intB_D ← inline(strrep(char(( dot(r(x, y, z) - rp(t),r(x, y, z)-rp(t)) )^(3/2)),"conjugate","")) 12. Hacer intB ← intB_N(t,x,y,z)/intB_D(t,x,y,z) 13. Hacer intBx ← dot(intB(t,x,y,z), ex) 14. Hacer intBy ← dot(intB(t,x,y,z), ey) 15. Hacer intBz ← dot(intB(t,x,y,z), ez)
16. Leer limi 17. Leer lims 18. Leer ene 19. Hacer A ← limi 20. Hacer B ← lims 21. Hacer N ← enero 22. Hacer H ← (B -A)/N 23. Hacer XI0 ← intBx(A) + intBx(B) 24. Hacer XY0 ← intBy(A) + intBy(B) 25. Hacer XZ0 ← intBz(A) + intBy(B) 26. Hacer XI1 ← 0.0 27. Hacer XY1 ← 0.0 28. Hacer XZ1 ← 0.0 29. Hacer XI2 ← 0.0 30. Hacer XY2 ← 0.0 31. Hacer XZ2 ← 0.0 32. Hacer NN ← N – 1 33. Hacer I ← 1
34. Repetir con I desde 1 hasta NN Hacer X ← A +( I * H) 35. 36. Si I MOD 2 es igual a 0 Hacer XI2 ← XI2 + intBx(X) 37. Hacer XY2 ← XY2 + intBy(X) 38. Hacer XZ2 ← XZ2 + intBz(X) 39.
35
40. Sino Hacer XI1 ← XI1 + intBx(X) 41. Hacer XY1 ← XY1 + intBy(X) 42. Hacer XZ1 ← XZ1 + intBz(X) 43. 44. Fin Si Hacer I ← I +1 45. 46. Fin de ciclo paso 34 47. Hacer Bx ← (XI0 + 2.0 * XI2 + 4.0 * XI1) * H / 3.0 48. Hacer By ← (XY0 + 2.0 * XY2 + 4.0 * XY1) * H / 3.0 49. Hacer Bz ← (XZ0 + 2.0 * XZ2 + 4.0 * XZ1) * H / 3.0
36
Solenoide { cx, cy, cz, I, limi, lims, ene son variables con punto flotante } 1. Hacer ex ← [1, 0, 0] 2. Hacer ey ← [0, 1, 0] 3. Hacer ez ← [0, 0, 1] 4. Hacer r ← [x, y, z] 5. Hacer rp ← [cos(t), sin(t), 0.05*t] 6. Hacer x ← cx 7. Hacer y ← cy 8. Hacer z ← cz 9. Hacer Km ← 1/(4*pi) 10. Hacer intB_N ← cross( diff( rp(t), t), r(x, y, z) - rp(t) ) 11. Hacer intB_D ← inline(strrep(char(( dot(r(x, y, z) - rp(t),r(x, y, z)-rp(t)) )^(3/2)),"conjugate","")) 12. Hacer intB ← intB_N(t,x,y,z)/intB_D(t,x,y,z) 13. Hacer intBx ← dot(intB(t,x,y,z), ex) 14. Hacer intBy ← dot(intB(t,x,y,z), ey) 15. Hacer intBz ← dot(intB(t,x,y,z), ez)
16. Leer limi 17. Leer lims 18. Leer ene 19. Hacer A ← lim i 20. Hacer B ← lims 21. Hacer N ← enero 22. Hacer H ← (B -A)/N 23. Hacer XI0 ← intBx(A) + intBx(B) 24. Hacer XY0 ← intBy(A) + intBy(B) 25. Hacer XZ0 ← intBz(A) + intBy(B) 26. Hacer XI1 ← 0.0 27. Hacer XY1 ← 0.0 28. Hacer XZ1 ← 0.0 29. Hacer XI2 ← 0.0 30. Hacer XY2 ← 0.0 31. Hacer XZ2 ← 0.0 32. Hacer NN ← N – 1 33. Hacer I ← 1
34. Repetir con I desde 1 hasta NN Hacer X ← A +( I * H) 35. 36. Si I MOD 2 es igual a 0 Hacer XI2 ← XI2 + intBx(X) 37. Hacer XY2 ← XY2 + intBy(X) 38. Hacer XZ2 ← XZ2 + intBz(X) 39. 40. Sino
37
Hacer XI1 ← XI1 + intBx(X) 41. Hacer XY1 ← XY1 + intBy(X) 42. Hacer XZ1 ← XZ1 + intBz(X) 43. 44. Fin Si Hacer I ← I +1 45. 46. Fin de ciclo paso 34 47. Hacer Bx ← (XI0 + 2.0 * XI2 + 4.0 * XI1) * H / 3.0 48. Hacer By ← (XY0 + 2.0 * XY2 + 4.0 * XY1) * H / 3.0 49. Hacer Bz ← (XZ0 + 2.0 * XZ2 + 4.0 * XZ1) * H / 3.0
38
Cylinder { ra, za, la, aa son variables con punto flotante } { er1, ep1, v son funciones } 1. Hacer er ← [1, 0, 0] 2. Hacer ep ← [0, 1, 0] 3. Hacer er1 ← [cos(p1), sin(p1), 0]; 4. Hacer ep1 ← [ -1*sin(p1), cos(p1), 0]; 5. Hacer ex ← [1, 0, 0] 6. Hacer ey ← [0, 1, 0] 7. Hacer ez ← [0, 0, 1] 8. Hacer v ← a*ep1(p1) 9. Hacer intn ← a*cross(v(z1, p1, a), (r*er -a*er1(p1)+(z-z1)*ez)) 10. Hacer intd ← inline(strrep(char(dot((r*er-a*er1(p1)+(z-z1)*ez),(r*er-a*er1(1)+(z-
z1)*ez)).^(3/2)), "conjugate", "")) 11. Hacer int ← (intn(a, p1, r, z, z1))/(intd(a, p1, r, z, z1)) 12. Hacer intx ← inline(strrep(char(dot(int(a, p1, r, z, z1),ex)), "conjugate", "")) 13. Hacer inty ← inline(strrep(char(dot(int(a, p1, r, z, z1),ey)), "conjugate", "")) 14. Hacer intz ← inline(strrep(char(dot(int(a, p1, r, z, z1),ez)), "conjugate", "")) 15. Hacer r ← ra 16. Hacer z ← za 17. Hacer L ← la 18. Hacer a ← aa 19. Hacer fx ← inline (strrep(strrep(strrep(char(eval("intx(a, p1, r, z, z1)")),"**",".^"), "*",".*"),
"/","./")) 20. Hacer Bx1 ← dblquad(fx,0,2*pi, -L/2,0) 21. Hacer Bx2 ← dblquad(fx,0,2*pi,0,L/2) 22. Hacer Bx ← Bx1+Bx2 23. Hacer fy ← inline(strrep(strrep(strrep(char(eval("inty (a, p1, r, z, z1)")),"**",".^"), "*",".*"),
"/","./")) 24. Hacer By1 ← dblquad(fy,0,2*pi, -L/2,0) 25. Hacer By2 ← dblquad(fy,0,2*pi,0,L/2) 26. Hacer By ← By1+By2 27. Hacer fz ← inline(strrep(strrep(strrep(char(eval("intz(a, p1, r, z, z1)")),"**",".^"), "*", ".*"),
"/","./")) 28. Hacer Bz1 ← dblquad(fz,0,2*pi, -L/2,0) 29. Hacer Bz2 ← dblquad(fz,0,2*pi,0,L/2) 30. Hacer Bz ← Bz1+Bz2
39
Spherical { rho, teta, ra son variables de punto flotante } 1. Hacer er ← [sin(t), 0, cos(t)] 2. Hacer et ← [cos(t), 0, -1*sin(t)] 3. Hacer ep ← [0, 1, 0] 4. Hacer er1 ← [sin(t1)*cos(p1), sin(t1)*sin(p1), cos(t1)] 5. Hacer et1 ← [cos(t1)*cos(p1), cos(t1)*sin(p1), -1*sin(t1)] 6. Hacer ep1 ←[ -1*sin(p1), cos(p1), 0] 7. Hacer ex ← [1, 0, 0] 8. Hacer ey ← [0, 1, 0] 9. Hacer ez ← [0, 0, 1] 10. Hacer v ← a*sin(t1)*ep1(p1) 11. Hacer intn ← (a.^2*sin(t1)*cross(v(t1, p1, a), (r*er(t) -a*er1(t1,p1)))) 12. Hacer intd ← inline(strrep(char(dot((r*er(t) -a*er1(t1, p1)),(r*er(t)-a*er1(t1, p1))).^(3/2)),
"conjugate", "")) 13. Hacer int ←intn(a, p1, r, t, t1)/intd(a, p1, r, t, t1) 14. Hacer intx ← inline(strrep(char(dot(int(a, p1, r, t, t1),ex)), "conjugate", "")) 15. Hacer inty ← inline(strrep(char(dot(int(a, p1, r, t, t1),ey)), "conjugate", "")) 16. Hacer intz ← inline(strrep (char(dot(int(a, p1, r, t, t1),ez)), "conjugate", "")); 17. Hacer a ← ra 18. Hacer r ← rho 19. Hacer t ← teta 20. Hacer fx ← inline(strrep(strrep(strrep(char(eval("intx(a, p1, r, t, t1)")),"**",".^"), "*",".*"),
"/","./")) 21. Hacer fy ← inline(strrep(st rrep(strrep(char(eval("inty(a, p1, r, t, t1)")),"**",".^"), "*",".*"),
"/","./")) 22. Hacer fz ← inline(strrep(strrep(strrep(char(eval("intz(a, p1, r, t, t1)")),"**",".^"), "*",".*"),
"/","./")) 23. Hacer Bx1 ← dblquad(fx,0,0.5,0,pi) 24. Hacer Bx2 ← dblqu ad(fx,0.5,(pi/2),0,pi) 25. Hacer Bx3 ← dblquad(fx,(pi/2),pi,0,pi) 26. Hacer Bx4 ← dblquad(fx,pi,((3*pi)/2),0,pi) 27. Hacer Bx5 ← dblquad(fx,((3*pi)/2),(2*pi),0,pi) 28. Hacer Bx ← Bx1+Bx2+Bx3+Bx4+Bx5 29. Hacer By1 ← dblquad(fy,0,0.5,0,pi) 30. Hacer By2 ← dblquad(fy,0.5,(pi/2),0,pi) 31. Hacer By3 ← dblquad(fy,(pi/2),pi,0,pi) 32. Hacer By4 ← dblquad(fy,pi,((3*pi)/2),0,pi) 33. Hacer By5 ← dblquad(fy,((3*pi)/2),(2*pi),0,pi) 34. Hacer By ← By1+By2+By3+By4+By5
35. Hacer Bz ← dblquad(fz,0,2*pi,0,pi)
40
EVIDENCIA DE FUNCIONAMIENTO Circular Loop
Current With General Shape
41
Solenoide
Cylinder
Spherycal
42
CONCLUSIONES El análisis numérico es una herramienta importante para resolver problemas cuya solución algebraica es complicada o imposible, y se pueden obtener resultados con buenas aproximaciones. Sin embargo como son soluciones aproximadas debemos de elegir un algoritmo que tenga poco error y sea rápido, en ese sentido no solo se trata de aplicar cualquier algoritmo, sino elegir entre los que existen la mejor opción. El análisis numérico se puede aplicar a muchas áreas, y es importante en el trabajo de científicos e ingenieros, en este trabajo se aplicó en el área de magnetismo, el cual tiene mucha importancia en nuestra vida diaria, por ejemplo en los discos duros de las computadoras ocupan el campo magnético para almacenar datos. También en los televisores que ocupaban un tubo de rayos catódicos, aunque ha quedado en desuso esa tecnología pero fue una aplicación importante durante muchos años. Sin embargo dependiendo de la aplicación del campo magnético que se estudie podemos obtener integrales bastante complejas que se pueden resolver fácilmente y con cierta rapidez a través de un método numérico. Las integrales simples se resolvían en cuestión de segundos, sólo se tuvo un problema en las integrales dobles, que eran un poco más lentas por la cantidad de operaciones que se debía realizar, en la documentación oficial encontramos lo siguiente al respecto de la lentitud de ese algoritmo: “Dicho enfoque funciona, pero es bastante lento, y ese problema aumenta
exponencialmente con la dimensionalidad de la integral. Otra solución posible es utilizar la Colocación Ortogonal como se describe en la sección anterior (ver Colocación Ortogonal).” Con colocación ortogonal era más rápido pero no lo ocupamos pues no se
hacía uso de ninguno de los métodos numéricos (vistos en clases o investigados), entonces se utilizó la funcion dblquad, que hace uso de la cuadratura gaussiana, uno de los métodos más exactos a la hora de realizar integración numérica. Como conclusión final queda agregar que es importante contribuir al desarrollo del software libre, ya que GNU/Octave es una herramienta maravillosa, pero comparado ante Matlab fue mucho más lento desarrollar la interfaz gráfica, por supuesto que es un software hecho por una empresa con millones de dólares y cuyas licencias profesionales llegan a valer más de $2000.00, por otro lado el software libre se mantiene a través de donaciones hechas por otros usuarios y los aportes que cada uno hace a la comunidad, tanto de código, documentación, testeo, entre otras. Sin embargo ha tenido bastante avance, antes ni siquiera contaba con una intefaz gráfica, solo era un programa de consola, cosa con la que ahora si se cuenta, y es un clon bastante bueno de uno de los programas más utilizados por cientificos e ingenieros como lo es Matlab.
43
RECOMENDACIONES Cuando se desee integrar con múltiples variables usando una de las funciones que provee Octave, tal es el caso de "dblquad”, se recomienda dividir los límites de
integración de un integrando con valores cercanos entre sí para un correcto funcionamiento del programa, a su vez manteniendo constante los límites del otro integrando.
44
BIBLIOGRAFÍA [1] Oliver Heimlich.(2016-2017).Numerical Integration.(GNU Texinfo 6.5).[Online].Disponible: https://octave.org/doc/v4.2.2/Numerical-Integration.html . [2] Young Hugh D. and Roger A. Freedman. “Fuentes de Campo Magnético” en Física
Universitaria con física moderna volumen 2: Pearson Education, México, 2009. [3] Hassani, Sadri. “Integration in Magnetism”, in Mathematical Methods Using Mathematica:
For Students of Physics and Related Fields, New York: Springer-Verlag, 2000,pp.110-123. [4] Álvarez Herrera, Javier.(26/08/2011).Integración Númerica.(156.35.0.0/16).[online].Disponible :https://www.unioviedo.es/compnum/laboratorios_ web/Laborat10_inte/Laborat07b_integracion.html .
45
GLOSARIO Algoritmo: conjunto ordenado de operaciones sistemáticas que permite hacer un cálculo y hallar la solución de un tipo de problemas. Campo magnético: espacio en que se hace sensible una fuerza determinada. Carga eléctrica: es una propiedad física intrínseca de algunas partículas subatómicas que se manifiesta mediante fuerzas de atracción y repulsión entre ellas a través de campos electromagnéticos. Coordenadas cilíndricas: las coordenadas cilíndricas son una extensión del sistema de coordenadas polares al espacio tridimensional. Corriente eléctrica: es el flujo de carga eléctrica que recorre un material. Cuadratura: en matemáticas, la cuadratura es un término histórico con el que se denomina la determinación del área de una figura. Dblquad: función para integrar con múltiples variables en Octave. Espira: un conductor cerrado plano se llama espira. Fórmula: es una secuencia o cadena de caracteres cuyos símbolos pertenecen a un lenguaje formal, de tal manera que la expresión cumple ciertas reglas de buena formación y que admite una interpretación consistente en alguna área de la matemática y en otros sistemas formales. Ésta tiene la finalidad de expresar una relación general entre los términos expresados en la fórmula. Función: en Octave y otros software una función es una entidad que recibe una serie de argumentos como entradas, opera sobre ellos y luego devuelve ciertos resultados. Integración numérica: en análisis numérico, la integración numérica constituye una amplia gama de algoritmos para calcular el valor numérico de una integral definida Integral: básicamente, una integral es una generalización de la suma de infinitos sumandos, infinitamente pequeños. Interfaz gráfica de usuario: La interfaz gráfica de usuario, conocida también como GUI (del inglés graphical user interface), es un programa informático que actúa de interfaz de usuario, utilizando un conjunto de imágenes y objetos gráficos para representar la información y acciones disponibles en la interfaz. Su principal uso, consiste en proporcionar un entorno visual sencillo para permitir la comunicación con el sistema operativo de una máquina o computador. Interpolación: en el subcampo matemático del análisis numérico, se denomina interpolación a la obtención de nuevos puntos partiendo del conocimiento de un conjunto discreto de puntos.
46
