ENSEÑANZA REVISTA MEXICANA DE FÍSICA E 53 (1) 56–66 JUNIO 2007

Animaciones en Matlab y maple de ecuaciones diferenciales parciales

de la fı́sica-matemática
G.M. Ortigoza Capetillo
Facultad de Matemáticas, Universidad Veracruzana, Zona Universitaria,
Apartado Postal 270, 91090 Xalapa, Ver.
e-mail: [email protected]
Recibido el 22 de mayo de 2006; aceptado el 23 de agosto de 2006
En este trabajo se presentan soluciones exactas de ecuaciones diferenciales parciales que dependen del tiempo; estas soluciones son de la
forma u(x, t), con x ∈ Rn , n = 1, 2, 3. Las gráficas de las soluciones a diferentes tiempos permiten la creación de animaciones de las
soluciones. Se muestra de manera general la forma de crear animaciones en Maple y Matlab. Estas animaciones pueden utilizarse como
herramienta didáctica para presentar fenómenos fı́sicos como son: la propagación de ondas de un medio a otro, superposición de ondas,
difusión, etc; ası́ mismo pueden usarse para despertar el interés de los estudiantes por el estudio de las ecuaciones diferenciales parciales y
sus aplicaciones. Para las animaciones se eligió un subconjunto importante de ecuaciones de la fı́sica matemática, entre las que se cuentan: la
ecuación del transporte, la ecuación de ondas (vibración de cuerdas y membranas, problema de transmisión), las ecuaciones de Klein Gordon,
Korteweg de Vries (no lineal), del calor y de Maxwell. Brevemente se describen algunas de las técnicas de solución analı́tica de edps como
son: escalamiento, método de caracterı́sticas, separación de variables, etc. Más aun, el contar con soluciones analı́ticas puede ser útil para la
verificación de implementaciones numéricas.

Descriptores: Enseñanza de la fı́sica; herramientas didácticas; ecuaciones diferenciales parciales.

In this work we present some exact solutions of time dependent partial differential equations (pdes); these solutions have the general form
u(x, t), with x ∈ Rn , n = 1, 2, 3. The plots of the solutions at different times allow us to create animations of the solutions. We show in a
general framework how to make animations in Maple and Matlab. These animations can be used as a didactic tool in order to introduce some
physical phenomena such as: wave propagation, superposition, transmission from one medium to another, diffusion, etc. They can also be
used to motive the students to the study of partial differential equations and its applications. A representative subset of differential equations
of mathematical physics was chosen that includes: the transport equation, wave equation, heat equation and equations of Klein Gordon,
Korteweg de Vries, and Maxwell. We briefly present some of the analytical methods for the solutions of pdes: scaling, characteristics and
separation of variables. Finally exact solutions can be very useful for code testing in numerical implementations.

Keywords: Physics education; education aids; partial differential equations.

PACS: 02.30Jr; 01.40Fk; 01.50Fr; 01.50Ht

1. Introducción de códigos para realizar animaciones en Matlab y Maple; pa-

Recientemente ha surgido un creciente interés por incorpo- ra que el lector, de acuerdo a sus propios intereses, pueda
rar, a la enseñanza de conceptos fı́sicos y matemáticos, el uso crear sus propias animaciones. El trabajo está organizado de
de nuevas tecnologı́as, como son: paquetes computacionales la siguiente manera: en la Sec. 2 se presentan los elementos
de cálculos simbólicos, animaciones y ambientes gráficos. básicos para realizar animaciones en 1, 2 y 3 dimensiones en
Ası́ podemos encontrar algunas fuentes bibliográficas don- Matlab y Maple; en la Sec. 3 introducimos las animaciones
de se estudia la fı́sica clásica , las ecuaciones diferenciales de ejemplos de una dimensión: la ecuación de transporte, la
y el análisis numérico usando los paquetes Maple, Matlab y ecuación de onda (sin fronteras, caso de intervalo finito, y un
Mathematica [3, 6–8, 11–13]. problema de transmisión), la ecuación de Klein-Gordon y se
concluye la sección con la ecuación de Korteweg de Vries, el
Estas herramientas computacionales han mostrado su uti- cual es un ejemplo de una ecuación no lineal. En la Sec. 3 se
lidad al presentar conceptos fı́sicos y matemáticos en forma muestran los ejemplos de dimensiones 2 y 3; entre los bidi-
amena y novedosa. Mediante el uso de animaciones se pue- mensionales se presentan el enfriamiento de una placa metáli-
de capturar la atención de los estudiantes y despertar su in- ca cuadrada y las vibraciones de una membrana circular, para
terés por incursionar en áreas de la fı́sica y las matemáticas los ejemplos tridimensionales se muestran la propagación de
aplicadas, tales como: mecánica celeste, teoria de vibracio- ondas electromagnéticas generadas por un dipolo eléctrico y
nes, mecánica de fluidos, ecuaciones diferenciales ordinarias, el enfriamiento de un cubo metálico. Finalmente incluimos
ecuaciones diferenciales parciales (edps), métodos numéri- algunas conclusiones.
cos, etc.
En este trabajo se presentan animaciones en Matlab y Ma- 2. Ejemplos de códigos para crear animacio-
ple de soluciones exactas de edps que dependen del tiempo, nes en Matlab y Maple
es decir, funciones de la forma u(x, tk ) con x ∈ Rn , tk ∈
R>0 ; se introducen de manera sencilla algunas de las técni- En esta sección mostraremos el uso de algunos comandos
cas más comunes de solución de edps y se describen ejemplos básicos para crear animaciones en Matlab y Maple.
 ANIMACIONES EN MATLAB Y MAPLE DE ECUACIONES DIFERENCIALES PARCIALES DE LA FÍSICA-MATEMÁTICA 57

2.1. Código para animaciones unidimensionales x=-2:0.1:10; nos define los valores de x. La ani-
mación se obtiene como una secuencia de los valo-
Supongamos que deseamos hacer la animación de la fun- res de u a diferentes tiempos, todos ellos graficados
2
ción u(x, t) = e−(x−2t) , la cual corresponde a la solución en un mismo marco; para ello utilizamos el comando
del problema de valor inicial: set(gca,’nextplot’,’replacechildren’); el cual mantiene las
2
ut + 2ux = 0, u(x, 0) = e−x . propiedades de los ejes. Para crear nuestra animación (se-
cuencia de gráficos) usamos un ciclo for dentro del cual
Veamos cómo se realiza esta animación en Matlab [16]: graficamos la solución al tiempo t(i) mediante el comando
La función de dos variables u(x, t) puede animarse grafi- plot(u,x) y guardamos el cuadro en una matriz M con la ayu-
cando en un mismo marco (plano x-u) sus valores a diferentes da del comando M(i)=getframe; en el transcurso de la eva-
tiempos. Ası́ definimos t=0 :0.05:6 el vector de los diferen- lución de este ciclo for, las gráficas a diferentes tiempos han
tes tiempos que toma como tiempo inicial 0, como tiempo sido guardadas en la matriz M y hemos visualizado la anima-
final 6 e incremento en tiempo de tamaño 0.05. Graficaremos ción. El siguiente código (mfile) muestra la implementación
en espacio (digamos en el intervalo [-2,10]), ası́ el vector de los comando anteriormente descritos:

% Animacion de la solucion de la ecuacion del transporte

% U_t + 2 U_x = 0
% con condicion inicial U(x,0)=exp(-xˆ2)
x=-2:0.1:10; t=0:0.05:6; set(gca,’nextplot’,’replacechildren’);
axis([-2 10 0 1]); for i=1:length(t)
u=exp(-(x-2*t(i)).ˆ2);
plot(x,u);
M(i)=getframe;
end
Para ver nuestra animación, abrimos el editor de Matlab y copiamos las instrucciones anteriores, grabamos el archivo
con extensión m, digamos ejemplo.m; para la ejecucción introducimos ejemplo en la ventana de comandos de Matlab. Una vez
ejecutado este código veremos nuestra animación, y si usamos el comando movie(M,0) ejecutaremos una vez más la animación.
Realizar la misma animación en Maple [14, 15] resulta más sencillo, para ello usamos las instrucciones. with(plots);
animate(plot, [ exp(-(x-2*t)ˆ2), x = -2 ..10], t= 0 .. 6);
Tecleando estas instrucciones en la ventana de comandos de Maple se generará una gráfica, damos un click en el botón
derecho del del ratón sobre ella y elegimos del menú Animate seguido por Play y de esta manera se obtendrá la animación.
La instrucción with(plots) le dice a Maple que debe usar el paquete plots, note que la forma de definir el espacio está dada
por x = −2..10 mientras que los valores para el parámetro t se definen como t = 0..6. Aquı́ el incremento del tiempo no fue
definido y Maple toma un valor por default, el cual hará nuestra animación más rápida o lenta (dependiendo de la velocidad de
procesamiento de nuestro sistema).

2.2. Código para animaciones bidimensionales

Mostraremos ahora un ejemplo de un código para simulaciones en 2 dimensiones; para ello consideremos la función
2
u(x, y, t) = e−π t sin(2πx) sin(2πy) la cual resuelve el problema de valor inicial con condiciones en la frontera:
µ ¶
1 ∂2u ∂2u
ut = + 2 (1)
8 ∂x2 ∂y
u(0, y, t) = u(1, y, t) = 0 u(x, 0, t) = u(x, 1, t) = 0 t>0
u(x, y, 0) = sin(2πx) sin(2πy). (2)

La animación en Matlab se consigue mediante un código muy similiar al caso unidimensional, a continuación mostramos un
código bidimensional donde mediante el uso del simbolo %, se han incluido comentarios breves acerca de los comandos más
importantes.
% Animacion de la solucion de la ecuacion del calor 2d
% U_t = 1/8(U_xx + U_yy)
% en el dominio 0 < x, y < 1
% condicion inicial u(x,y,0)=Sin(pix) Sin(piy)
% para 0<t<3

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66

58 G.M. ORTIGOZA CAPETILLO

[x y]=meshgrid(0:.01:1,0:.01:1);
%define una malla para [0,1]x[0,1] con incremento en x y y de 0.01
title(’Animacion de la temperatura’);
set(gca,’nextplot’,’replacechildren’); caxis manual;
% permite que todos los gráficos usen los mismos lı́mites en colores.
caxis([-1 1]);
% define los valores maximos para los limites como -1 y 1
axis equal;
%usa la misma escala para los ejes x y
t=0:0.05:3;
%vector de valores para los diferentes tiempos
for j=1:length(t) z=exp(-1/8*piˆ2.*t(j)).*sin(2*pi*x).*sin(2*pi*y);
%evaluacion de la funcion de dos variables
axis off;
%remueve los ejes
pcolor(x,y,z);
%grafica en dos dimensiones
shading interp;
%interpolacion de colores
colorbar;
% agrega la barra de colores
M(j) = getframe(gcf);
% captura los graficos y los guarda en la matriz M
end

Para obtener la misma animación en Maple, usamos la

siguiente instrucción:
realidad tenemos una gráfica en 3d pero la visualizamos
animate(plot3d, [exp(-Pi*Pi*t)*sin(2*Pi*x)*sin(2*Pi*y), x en el plano xy, esto se logra cambiando la orientación
= 0 .. 1, y = 0 .. 1], t = 0 .. .0.2, orientation = [-90, 0], style = orientation[θ, φ,] , patchnogrid usa un reticulado para colo-
patchnogrid, axes = boxed, shading = zhue); rear sin mostrarlo, axes muestra los ejes coordenados y final-
Expliquemos un poco esta instrucción; el comando básico pa- mente shading=zhue ajusta los colores de acuerdo al valor
ra hacer animaciones en Maple es animate y su sintaxis es la de la función en cada punto del reticulado. Una vez ejecutada
siguiente: esta instrucción la figura obtenida se anima con el ratón de la
animate(plotcommand, plotargs, t=a..b, options); misma manera que para el caso unidimensional.
aquı́ plotcommand es el procedimiento Maple con el cual se
2.3. Código para Animaciones Tridimensionales
obtiene un gráfico 2D o 3D; plotargs es la lista de argumen-
tos del comando de gráfico, t es el parámetro mediante el cual Para presentar los ejemplos de códigos de animaciones tridi-
se consigue la animación el cual varı́a entre los valores a y b mensionales usaremos la función:
y finalmente options especifica opciones del comando para
2
graficar o de animación. U (x, y, z, t) = e−π t sin(πx) sin(πy) sin(πz),
Para nuestra animación el comando para graficar es
que corresponde a la solución exacta del problema de la
plot3d, sus argumentos aparecen entre corchetes cuadrados
Sec. 4.4. El parámetro t corresponde al tiempo, el dominio
la función de dos variables espaciales y del parámetro tiempo
espacial es Ω = [0, 1]3 y para la visualización usaremos cor-
(se incluye el rango de definición de las variables espaciales);
tes definidos por la intersección de Ω con planos paralelos. El
con respecto a las opciones incluidas mencionamos que en
siguiente código Matlab es muy similar a los dados anterior-
mente.

[x,y,z] = meshgrid(0:.05:1,0:0.05:1,0:.05:1);
set(gca,’nextplot’,’replacechildren’); caxis manual; caxis([0,1]);
t=0:0.01:0.5; for j=1:length(t)
v=exp(-pi*pi*t(j)).*sin(pi*x).*sin(pi*y).*sin(pi*z);
h = slice(x,y,z,v,[0.25 .5 .75],[],[0]);
caxis([0,1]); alpha(’color’)
set(h,’EdgeColor’,’none’,’FaceColor’,’interp’,...

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66

 ANIMACIONES EN MATLAB Y MAPLE DE ECUACIONES DIFERENCIALES PARCIALES DE LA FÍSICA-MATEMÁTICA 59

’FaceAlpha’,’interp’)
alphamap(’rampdown’) alphamap(’increase’,.1) colormap(hsv);
colorbar; F(j)=getframe(gcf); end;
Para la visualización los cortes se obtienen mediante la instrucción h = slice(x,y,z,v,[0.25 .5 .75],[ ],[0]) que corresponde
a Ω intersectado con x = 0.25, x =0.5 y x =0.75 ( [16]) . La instrucción set se usa para definir las opciones del gráfico h, alpha
y alphamap para definir las propiedades de transparencia del gráfico mientras que colormap define los colores empleados.
Para la misma animación en Maple usamos el siguiente código:
with(plots):N := 200:anim := array(1 .. N): for i to 25 do t := i/N:
P1[i] := plot3d(.25, x = 0 .. 1, y = 0 .. 1, color =
exp(-Pi*Pi*t)*sin(Pi*x)*sin(Pi*y)*sin(.25*Pi), style = patchnogrid);
P2[i] := plot3d(.5, x = 0 .. 1, y = 0 .. 1, color =
exp(-Pi*Pi*t)*sin(Pi*x)*sin(Pi*y)*sin(.5*Pi), style = patchnogrid);
P3[i] := plot3d(.75, x = 0 .. 1, y = 0 .. 1, color =
exp(-Pi*Pi*t)*sin(Pi*x)*sin(Pi*y)*sin(.75*Pi), style = patchnogrid);
anim[i] := display([P1[i], P2[i], P3[i]]) od;
display([seq(anim[i], i = 1 .. 25)], insequence = true, axes =
boxed, transparency = .5)

En este caso hemos utilizado un arreglo para guardar la

solución a diferentes tiempos; la solución se compone de Ası́,
los gráficos P 1, P 2 y P 3 que corresponden a los diferen-
tes cortes mediante los planos z =0.25, z =0.5 y z =0.75. v(t) = v(0) v(t) = u(x(0), 0) = u(x0 , 0) = f (x0 )
El comando display nos da como resultado una gráfica, la Por lo tanto la solución del problema de valor inicial (3) es
cual se puede animar con los botones del ratón. Se ha utili-
zado el comando color para asignar el color correspondiente u(x, t) = f (x − ct).
a cada punto del espacio tridimensional y nótese que para
nuestra animación el incremento del tiempo está dado por La solución se obtiene desplazando ct unidades a la dere-
t = i/N [14, 15]. cha la condición inicial (c>0). En la animación considera-
2
mos c = 2 y f (x) = e−x . La Fig. 1 muestra la solución para
t = 0, 1, 2, 3. Observamos ası́ que la ecuación de transporte
3. Animaciones de edps en una dimensión tiene el efecto de mover la condición inicial a la derecha sin
distorsionarla.
3.1. La ecuación de transporte

Considere el problema de valor inicial 3.2. La ecuación de ondas (sin fronteras)

ut + cux = 0, u(x, 0) = f (x), c constante. (3) Considere las ecuaciones

Para resolver este problema utilizaremos el método de carac- utt = c2 uxx , u(x, 0) = f (x),
ut (x, 0) = g(x).
terı́sticas: (5)
La Ec. (3) puede escribirse como Usaremos en este ejemplo el método de cambio de varia-
∂u bles [2]. Definimos
(c, 1) · (ux , ut ) = 0 ~ =0
= ∇u · V (4)
~
∂V ξ = x − ct, η = x + ct;
Esta ecuación nos dice que u es constante en el plano x − t
~ = (c, 1). A estas con este cambio de variables las Ecs. (5) se transforman en
a lo largo de las rectas paralelas al vector V
curvas se les llama curvas Caracterı́sticas. uξη = 0
Considere la recta que al tiempo t = 0 corta al eje x en
x0 , su parametrización está dada por cuya solución esta dada por:
½
x(t) = x0 + c t u(ξ, η) = F (ξ) + G(η).
t = t.
Por lo tanto, u expresada en las variables x, t tiene la forma
La función
v(t) = u(x(t), t) u(x, t) = F (x − ct) + G(x + ct).
es constante; en efecto De acuerdo con las condiciones iniciales,
dv ∂u dx ∂u ∂u ∂u
= + = c+ =0 f (x) = F (x) + G(x), g(x) = −cF 0 (x) + cG0 (x)
dt ∂x dt ∂t ∂x ∂t

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66

60 G.M. ORTIGOZA CAPETILLO

En este caso tenemos que la condición inicial se divide en

dos ondas propagandose a la misma velocidad en direcciones
opuestas.

3.3. La ecuación de ondas en un intervalo finito (vibra-

ción de una cuerda)

Considere el problema de valor inicial con condiciones de

frontera
utt = a2 uxx , 0 < x < L, t > 0, (6)
u(0, t) = 0, u(L, t) = 0, 0 < t < ∞, (7)
u(x, 0) = f (x), ut (x, 0) = g(x), 0 ≤ x ≤ L. (8)
Para obtener la solución usaremos el método de sepación de
variables
Supongamos que u tiene la forma
F IGURA 1.
u(x, t) = X(x)T (t), (9)
al substituir (9) en (6) obtenemos
T 00 X 00
XT 00 = a2 X 00 T 2
= = λ, λ constante
a T X
½
T 00 − λa2 T = 0
(10)
X 00 − λX = 0.
Supongamos que λ < 0, λ = −β 2 . Las soluciones de (10)
están dadas por

T (t) = A sin(aβt) + B cos(aβt), (11)

X(x) = C sin(βx) + D cos(βx).
Usando las condiciones de frontera (7) tenemos que
X(0) = 0, X(L) = 0.
F IGURA 2. Ası́, X(0) = 0 implica que D = 0 y X(L) = 0 implica que
βL = nπ con n ∈ Z, con lo cual
y resolviendo para F y G obtenemos
³ nπ ´ ³ nπ ´
Zx Zx Tn (t) = A sin a t + B cos a t , (12)
f (x) 1 f (x) 1 L L
F (x)= − g(s)ds, G(x)= + g(s)ds. ³ nπ ´
2 2c 2 2c
0 0 Xn (x) = C sin x ,
L
Por lo tanto, la solución del problema de valor inicial (5) tiene y por superposición tenemos
la forma X∞ ³ nπ ´
Z
x+ct u(x, t) = sin x
1 1 n=1
L
u(x, t) = [f (x − ct) + f (x + ct)] + g(s)ds. h ³ nπ ´ ³ nπ ´i
2 2c
x−ct × An sin a t + Bn cos a t .
L L
A la expresión anterior se le conoce como fórmula Usando ahora las condiciones iniciales (8) tenemos que
de D’Alembert. [2] En la animación consideramos c=1, X∞ ³ nπ ´
2
f (x) = e−2(x−2) y g = 0, ası́ la solución es u(x, 0) = f (x) = Bn sin x ,
L
1 h −2(x−t−2)2 i n=1
2
u(x, t) = e + e−2(x+t−2) X∞
nπa ³ nπ ´
2 ut (x, 0) = g(x) = An sin x ,
L L
La Fig. 2 muestra la solución para t = 0, 1, 2, 3. n=1

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66

 ANIMACIONES EN MATLAB Y MAPLE DE ECUACIONES DIFERENCIALES PARCIALES DE LA FÍSICA-MATEMÁTICA 61

F IGURA 3.
F IGURA 4.
donde
ZL ³ nπ ´ Las condiciones para unir las soluciones son la continuidad
2
An = g(x) sin x dx, de u y ux en x = 0,
nπa L
0

ZL ³ nπ ´  H(cl t) = F (−cr t) + G(cr t), ∀t,
2 (14)
Bn = f (x) sin x dx. 
L L H 0 (cl t) = F 0 (−cr t) + G0 (cr t) ∀t.
0
Para la animación tomamos L = 2, a = 1, g ≡ 0 y
³ πx ´ 1 µ ¶ µ ¶ resolviendo estas ecuaciones para F y para H obtenemos
3πx 1 5πx
f (x) = sin + sin + sin .
2 2 2 4 2
cl − cr
Por lo tanto, la solución está dada por F (z) = G(−z),
³ πx ´ µ ¶ µ ¶ µ ¶ cl + cr
πt 1 3πx 3πt µ ¶
u(x, t)= sin cos + sin cos 2cl cr z
2 2 2 2 2 H(z) = G , z = cl t.
µ ¶ µ ¶ cl + cr cl
1 5πx 5πt
+ sin cos .
4 2 2
Por lo tanto la solución u se escribe como
La Fig. 3 muestra la solución para t = 0, 0.25, 0.5, 0.75. No-
tamos que los extremos de la cuerda están fijos y cualquier 
otro punto de la cuerda oscila verticalmente.  G(x + cr t) + RG(cr t − x), x > 0,
u(x, t) = (15)
 T G( cr (x + c t)), x ≤ 0.
3.4. Un problema de transmisión cl l

Considere dos resortes de diferentes materiales unidos en

x = 0, sujetos a una misma tensión T0 pero de diferentes Con R = (cl − cr )/(cl + cr ) y T = 2cl /(cl + cr ).
densidades ρl y ρr [3]. Ası́ las velocidades de propagación Para la animación consideramos cl = 1 cr = 3 y
son p p
cl = T0 /ρl , cr = T0 /ρr . 2
(x − 5)(x − 6)e−(x−5)
Supongamos que tenemos una onda incidente G(x+cr t) mo- G(x) = .
0.7635
viéndose a la izquierda con G(x) = 0, x ≤ 0. Además de
la onda reflejada F (x − cr t) obtendremos una onda trans-
mitida H(x + cl t) moviéndose a la izquierda. La solución La Fig. 4 muestra la solución para t = 0, 1, 2, 3. Inicialmente
estará dada por una onda que se mueve a la izquierda, al alcanzar el punto
 x = 0, se divide en una onda transmitida (hacia la izquierda)
 F (x − cr t) + G(x + cr t), x > 0 y una onda reflejada (hacia la derecha). Este ejemplo permite
u(x, t) = (13) mostrar el fenómeno de transmisión de ondas de un medio a

H(x + cl t), x ≤ 0. otro.

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66

62 G.M. ORTIGOZA CAPETILLO

3.6. Korteweg de Vries

Consideremos la ecuación no lineal

ut + 6uux + uxxx = 0. (19)

Buscamos soluciones de la forma u(x, t) = v(x−σt) = v(s)

(onda viajera [1]). Al substituir en (19) obtenemos

−σv 0 + 6vv 0 + v 000 = 0,

−σv + 3v 2 + v 00 = c1 ,
σv 2 v 02
− + v3 + = c1 v + c2 .
2 2

Si requerimos que v, v 0 , v 00 → 0 cuando t → ±∞, entonces

F IGURA 5. Pie de figura. c1 = 0 y c2 = 0. Resolviendo para v 0 en la última expresión
obtenemos
3.5. Ecuación de Klein-Gordon √
v 0 = ±v σ − 2v,
Consideremos la ecuación de Klein-Gordon
tomando el signo menos e integrando obtenemos
utt − uxx + m2 u = 0, m constant. (16)
v(s)
Z
Definamos el cambio de variables dz
s=− √ + c. (20)
ξ = t − x, η = t + x, z σ − 2z
0

y consideremos
U (ξ, η) = u(x, t). Considerando el cambio de variable z = (σ/2)sech2 (θ) te-
nemos que dz/dθ = −σsech2 (θ) tanh(θ) y después de inte-
Si u satisface (16), entonces U (ξ, η) satisface grar (20) obtenemos
m2
Uξη + U = 0. (17) 2 σ
4 s = √ θ + c, sech2 (θ) = v(s).
σ 2
Esta ecuación es invariante bajo la tranformación ξ → αξ,
η → α−1 η (método de escalamiento), ası́ proponemos una combinando estas expresiones tenemos que
solución de la forma
µ√ ¶
U (ξ, η) = f (z), z = ξη, σ σ
v(s) = sech2 (s − c) ,
2 2
teniendo ası́ que f (z) satisface
m2 y ası́ finalmente tenemos que la solución u(x, t) tiene la for-
zf 00 (z) + f 0 (z) + f =0 ma
4 µ√ ¶
√ σ 2 σ
y tomando s = m z obtenemos la ecuación de Bessel de u(x, t) = sech (x − σt − c) .
2 2
orden cero:
sf 00 + f 0 + sf = 0. La expresión anterior de la solución nos dice qué ondas via-
J0 (s) es una solución regular para s = 0. Ası́, proponemos jeras con velocidad de propagación σ tienen amplitud pro-
la solución porcional a σ. Ondas de mayor amplitud viajan con mayor
½ 1 √ velocidad. En la animación se consideran dos ondas con ve-
J0 (m( t2 − x2 )), x2 < t2 , locidades de propagación σ1 = 5 y σ2 = 2, la primera cen-
u(x, t) = 2 (18)
0, x2 > t2 . trada en x = 0 y la segunda en x = 5. La Fig. 6 muestra la
Para la animación consideramos m = 4. La Fig. 5 muestra solución para t = 0, 1, 2, 3, 4. Notamos que la onda de mayor
la solución para t = 0, 1, 2, 3, 4, 5. Esta solución se comporta amplitud (y mayor velocidad) alcanza y cruza a la de menor
de manera similar a la ecuación de ondas: dos ondas se pro- amplitud.
pagan en direcciones opuestas; sin embargo se tiene el efecto Las gráficas de las animaciones de esta sección fueron ob-
de fuente en el punto x = 0, donde se generan ondas que se tenidos modificando la función u(x, t), los vectores espacio
propagan a la derecha y a la izquierda. y tiempo x, t del código Matlab dado en la Sec. 2.1.

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66

 ANIMACIONES EN MATLAB Y MAPLE DE ECUACIONES DIFERENCIALES PARCIALES DE LA FÍSICA-MATEMÁTICA 63

substituyendo esta expresión en 21 tenemos

µ ¶
1 1
RΘT 00 = c2 R00 ΘT + R0 ΘT + 2 RΘ00 T , (24)
r r
es decir,
µ ¶
T 00 R00 1 R0 1 Θ00
= + + 2 = −λ2 . (25)
c2 T R r R r Θ
Para separar R y Θ consideramos:
r2 R00 rR0 Θ00
+ + + λ2 r2 = 0, (26)
R R Θ
r2 R00 rR0 Θ00
+ + λ2 r 2 = − = n2 . (27)
R R Θ
De las Ecs. (25) y (27) obtenemos
 2 00

 r R + rR0 + (λ2 r2 − n2 )R = 0,
F IGURA 6. 


Θ00 + n2 Θ = 0, (28)




 00
T + c2 λ2 T = 0,
cuyas soluciones están dadas por


 R(r) = a1 Yn (λr) + a2 Jn (λr),



Θ(θ) = b1 cos(nθ) + b2 sin(nθ), (29)





T (t) = c1 cos(λct) + c2 sin(λct).
Las condiciones iniciales (22) nos dan que

R(0) < ∞ ⇒ a1 = 0 ∧ R(1) = 0

⇒ λ = knm m = 1, 2, . . .

Ası́, por superposición, u está dada por

F IGURA 7. ∞ X
∞
X
u(r, θ, t) = Jn (knm r) cos(nθ)
n=0 m=1
4. Ejemplos en dos y tres dimensiones
× [Anm sin(knm ct) + Bnm cos(knm ct)] . (30)
4.1. La ecuación de ondas en coordenadas polares:
vibración de una membrana circular Consideraremos el caso

Considere el problema de valor inicial con valores en la fron- u(r, θ, 0) = f (r), ut (r, θ, 0) = 0,
tera:
de esta manera (30) se reduce a
µ ¶
2 1 1 ∞
X
utt = c urr + ur + 2 uθθ , (21)
r r u(r, t) = Am J0 (k0m r) cos(k0m ct). (31)
m=1
u(1, θ, t) = 0, 0 < t < ∞, 0 < θ < 2π, (22)
Por la condición inicial
u(r, θ, 0) = f (r, θ), ut (r, θ, 0) = g(r, θ). (23) ∞
X
f (r) = Am J0 (k0m r),
Usaremos nuevamente el método de separación de varia-
m=1
bles [4].
Supongamos que u tiene la forma Z1
2
Am = 2 rf (r)J0 (k0m r)dr.
J1 (k0m )
u(r, θ, t) = R(r)Θ(θ)T (t), 0

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66

64 G.M. ORTIGOZA CAPETILLO

Para la animación asumiremos

1
c = 1, u(r, θ, 0) = J0 (2.4r) + J0 (8.65r),
2
ut (r, θ, 0) = 0.

Por lo tanto la solución u es

1
u(r, t) = J0 (2.4r) cos(2.4t) + J0 (8.65r) cos(8.65t).
2
La Fig. 7 muestra la solución para t=0, 0.25, 0.5, 0.75, 1 y
1.25. Notamos que la membrana oscila verticalmente.

4.2. La ecuación del calor

Considere el problema de valor inicial con condiciones de

F IGURA 8.
frontera:
µ 2 ¶ 4.3. Ecuaciones de Maxwell : dipole eléctrico
∂ u ∂2u
ut = k + 2 , (32)
∂x2 ∂y Las ecuaciones de Maxwell en espacio libre en ausencia de
u(0, y, t)=u(1, y, t)=0, u(x, 0, t)=u(x, 1, t)=0, t>0, (33) cargas están dadas por
∂H ~
u(x, y, 0)=f (x, y), 0 ≤ x ≤ 1, 0 ≤ y ≤ 1. (34) µ0 = −∇ × E,~ ∇·E ~ = 0,
∂t
Este problema de valor inicial modela la variación de la tem- ∂E ~
²0 ~
= ∇ × H, ∇·H ~ = 0. (36)
peratura de una placa cuadrada, cuya frontera se mantiene a ∂t
temperatura cero. Usando separación de variables [5] obte- En este ejemplo usaremos el método de potenciales [9]. Su-
nemos: pongamos que el dipolo es pequeño y que la corriente I es
uniforme sobre su longitud. Sea dl la longitud del dipolo.
u(x, y, t) Consideremos un sistema en coordenadas esféricas. Supon-
∞ X
X ∞ gamos que la corriente varı́a en forma sinusoidal y ası́ el vec-
2
+n2 )π 2 t
= Bmn e−k(m sin(mπx) sin(nπy), (35) tor potencial está dado por
m=1 n=1 Z
µ0 e−jβr
A(~r) = J(r) dv.
donde 4π r
V
∞ X
X ∞
Supongamos además que la corriente está encerrada en
f (x, y) = Bmn sin(mπx) sin(nπy),
un cable cuya sección transversal ds es pequeña. Ası́,
m=1 n=1
Jdv = Jdsdl = Idlẑ para un dipolo localizado en el origen
Z1 Z1 y orientado en la dirección del eje ẑ, el vector potencial se
Bmn = 4 f (x, y) sin(mπx) sin(nπy)dxdy. reduce a
µ0 Idl e−jβr
0 0 A(~r) = ẑ,
√ 4π r
Para la animación tomamos donde β = ω µ0 ²0 .
El vector A en coordenadas esféricas está dado por
1
k= f (x, y) = sin(2πx) sin(2πy). µ0 Idl e−jβr µ0 Idl e−jβr
8 A= cos(θ)r̂ − sin(θ)θ̂.
4π r 4π r
Ası́, la solución u está dada por Los vectores de fase magnético y eléctrico están dados por
2 1 1
u(x, y, t) = e−π t sin(2πx) sin(2πy). H= ∇ × A, E= ∇ × H,
µ0 jω²0
La Fig. 8 muestra la temperatura para t=0, 0.2, 0.4, 0.6, Los vectores magnético y eléctrico se obtienen como
0.8, 1. La temperatura toma el valor mı́nimo -1 y máximo 1. ~ = Re{ejωt H},
H ~ = Re{ejωt E},
E
Conforme pasa el tiempo los puntos más frı́os (azul) aumen- · ¸
tan su temperatura (rojo) y los más calientes disminuyen has- ~ Idl cos(ωt − βr) β sin(ωt − βr)
H= − sin(θ)φ̂,
ta alcanzar el estado estacionario de temperatura cero. 4π r2 r

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66

 ANIMACIONES EN MATLAB Y MAPLE DE ECUACIONES DIFERENCIALES PARCIALES DE LA FÍSICA-MATEMÁTICA 65

4.4. Una solución de la ecuación del calor en 3 dimen-

siones

Consideremos el problema de valor en la frontera con condi-

ciones iniciales
Ut = ν(Uxx + Uyy + Uzz ) en Ω,
U (x, y, z, t)|∂Ω = 0, U (x, y, z, 0) = g(x, y, z), (37)
donde asumiremos que Ω=[0, 1]3 . Usando separa-
ción de variables de manera similar al ejemplo bidi-
mensional, y suponiendo que la condición inicial es
g(x, y, z) = sin(πx) sin(πy) sin(πz), concluimos que
2
U (x, y, z, t) = e−π t sin(πx) sin(πy) sin(πz)
es la solución de (37) con ν = 1/3. La Fig. 10 muestra la tem-
peratura a tiempos: t=0, 0.5, 0.1 y 1.5. Conforme el tiempo
F IGURA 9.
avanza el cubo se enfria hasta alcanzar el estado estacionario
de temperatura cero. Para la visualización se usaron cortes
con 3 planos.
Los gráficos de las Secs. 4.1, 4.2 y 4.3 fueron obtenidos
con códigos Matlab muy similares al código bidimensional
dado en la Sec. 2.2, mientras que el ejemplo tridimensional de
la Sec. 4.4 fue obtenido con el código Matlab de la Sec. 2.3.

5. Conclusiones
En este trabajo se han presentado animaciones en Matlab de
soluciones exactas de edps de la forma u(x, tk ) con x ∈
Rn , tk ∈ R>0 . Se ha seleccionado un conjunto de las ecua-
ciones más representativas de la fı́sica matemática. Estas ani-
maciones pueden utilizarse para explicar fenómenos fı́sicos
como superposición de ondas, transmición de ondas de un
medio a otro, etc; ası́ como para despertar el interés de los es-
tudiantes y encaminarlos en el área de las edps y sus aplica-
F IGURA 10. ciones. Se han incluido códigos en Matlab y Maple para crear
animaciones en 1, 2 y 3 dimensiones, ası́ como brevemente
se describen algunas de las técnicas de solución de las edps:
· ¸
Idl β cos(ωt − βr) sin(ωt − βr) escalamiento, cambio de variables, método de caracterı́sticas,
~ =
E + cos(θ)r̂
2π²0 ω r2 r3 separación de variables, etc.
· La mayoria de las veces las edps halladas en las aplicacio-
Idl β sin(ωt − βr) β cos(ωt − βr) nes no tiene soluciones analı́ticas, en cuyo caso se tienen que
+ − +
4π²0 ω r r2 emplear métodos numéricos de aproximación. Contar con al-
¸ gunas soluciones exactas nos permite verificar la veracidad
sin(ωt − βr)
+ sin(θ)θ̂ de los métodos numéricos ( [10]). Por otra parte, el uso de
r3
animaciones de ecuaciones puede ser útil en el modelado:
Para la animación consideramos Idl = 4π, β = ω = 2π, modificando determinados parámetros en las ecuaciones o en
²0 = µ0 = 1 y por simetrı́a usamos la proyección sobre el sus soluciones pueden ajustarse las animaciones a una me-
eje xy. La Fig. 9 muestra la superposición de dos dipolos pa- jor representación de un fenómeno fı́sico, obteniendo de esta
ra t=0.5, 1.0, 1.5, y 2.0, donde se ha graficado la intensidad manera modelos matemáticos que describan con mayor exac-
del campo magnético. Observamos cómo se generan las on- titud un fenómeno.
das electromagnéticas esféricas propagándose hacia afuera y
se superponen de manera similar a las ondas en el agua. Por Agradecimientos
simetrı́a se graficó un corte por un plano, pero en realidad
éste es un ejemplo de una solución exacta de las ecuaciones Trabajo realizado con apoyo de proyecto Promep
de Maxwell en 3 dimensiones [10]. 103.5/05/1955.

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66

66 G.M. ORTIGOZA CAPETILLO

1. L.C. Evans, Partial Differential Equations (Graduate Studies in 9. D.J. Griffiths, Introduction to Electrodynamics (Prentice Hall,
Mathematics, American Mathematical Society, 1998). 1999).
2. J D. Logan, Applied Partial Differential Equations (Springer- 10. G. Ortigoza, The Runge-Kutta Discontinuous Galerkin method
Verlag, Undergradute Texts in Mathematics, 1998). for Maxwell equations, Ph. D. thesis, 2003.
3. J. Cooper, Introduction to Partial Differential Equations with
11. K.A. Lonngren and Sava Savov Fundamentals of Electromag-
Matlab (Birkhäuser, 1998).
netics with MATLAB (SciTech Publishing, Incorporated, 2005).
4. S.J. Farlow, Partial Differential Equations for Scientists and en-
gineers (Dover, 1993). 12. F. Wang, Physics with MAPLE: The Computer Algebra Resour-
5. P. Duchateau and D.W. Zachmann, Partial Differential Equa- ce for Mathematical Methods in Physics (Wiley, John and Sons
tions (Schaum´s Outline series, McGraw-Hill,1986). Incorporated, 2006).
6. I.P. Stavroulakis and S.A. Tersian, Partial Differential Equa- 13. Gerd Baumann , Mathematica for Theoretical Physics: Classi-
tions: An Introduction with Mathematica and Maple (World cal Mechanics and Nonlinear Dynamics (Springer-Verlag New
Scientific Publishing Company, 2004). York, 2005).
7. L. Elden, L. Wittmeyer-Koch, and H.B. Neilsen, Introduction to 14. Maple 10 ,user manual, maplesoft, 2005.
Numerical Computation - Analysis and MATLAB Illustrations
(Studentlitteratur, 2004). 15. J. Putz,Maple Animation (CRC Press, May 2003).
8. M.L. Abell and J.P. Braselton, Differential Equations with Mat- 16. D. Hanselman and B. Littlefield, Mastering Matlab 7, (Pear-
hematica (Elsevier Science & Technology Books, 2004). son/Prentice Hall, 2005).

Rev. Mex. Fı́s. E 53 (1) (2007) 56–66