Apunte Matematica D
Apunte Matematica D
Apunte Matematica D
Departamento de Matemática
Facultad de Ingenierı́a Quı́mica
Universidad Nacional del Litoral
16 de marzo de 2016
Capı́tulo 1
El error relativo mira el error como una proporción del valor exacto. En otras
palabras, es como si tuviera en cuenta las unidades de medida
En general, los métodos numéricos permiten hallar una aproximación p̂ de una solución
exacta p y proporcionan una estimación para el error absoluto Ep .
Consideremos el caso en que sabemos que el error absoluto Ep es menor o igual a 0.02.
Entonces
2
Es decir, una aproximación p̂ y una cota para el error absoluto nos permiten establecer
un intervalo donde está el valor exacto p.
Para el razonamiento que sigue usaremos las siguientes propiedades del valor absoluto:
|x + y| ≤ |x| + |y|, |x − y| ≤ |x| + |y|, |x| − |y| ≤ |x − y|, |y| − |x| ≤ |x − y|.
1 1
|p̂| − |p| ≤ |p̂ − p| = Ep ≤ ep =⇒ |p| ≥ |p̂| − ep =⇒ ≤ ,
|p| |p̂| − ep
y luego
|p − p̂| ep
Rp = ≤ .
|p| |p̂| − ep
ep
Cuando ep es muy pequeño en relación a |p̂| (un 10 % o menos), entonces es muy
|p̂| − ep
ep
aproximado a .
|p̂|
Muchos algoritmos como los que veremos más adelante se detienen cuando se cumple
ep
que es menor a una tolerancia preestablecida para el error relativo.
|p̂|
En general, un método numérico para calcular soluciones de ecuaciones provee suce-
sivas aproximaciones p̂ y estimaciones del error absoluto ep .
Si se desea una aproximación con un error absoluto menor a una tolerancia tol,
detendremos el procedimiento cuando se cumpla ep < tol.
Si se desea una aproximación con un error relativo menor a una tolerancia tolr,
ep
detendremos el procedimiento cuando se cumpla < tolr.
|p̂|
>> sqrt(6)
ans =
2.4495
>> 2.4495^2
ans =
6.0001
√
¿Qué ocurrió? Ocurre que 2.4495 no es igual a 6 sino sólo una aproximación. Es im-
portante notar que aunque MATLAB muestra que la respuesta a sqrt(6) es 2.4495, en
realidad calculó la raı́z de 6 con 16 dı́gitos exactos. Veamos:
>> sqrt(6)
ans =
3
2.4495
>> ans^2
ans =
6.0000
Para evitar escribir todo el número o copiar y pegar, usamos la variable ans que almacena
el último resultado calculado (que no es asignado a otra variable) con todos los dı́gitos
calculados. Por eso parece que el cálculo es exacto. Sin embargo los cuatro ceros en 6.0000
indican que el resultado no es exactamente 6. Veamos el error:
>> abs(ans-6)
ans =
8.8818e-16
>> abs(6-ans)
ans =
8.881784197001252e-16
>> x = 0.25
x =
0.250000000000000
>> 1+x-1
ans =
0.250000000000000
4
>> x = 2e-20
x =
2.000000000000000e-20
>> 1+x-1
ans =
0
¿Eh? ¿Qué ocurrió? Sucede que cuando la computadora suma 1 más 2 × 10−20 la suma
le da 1 porque redondea al número más cercano que puede representar con 16 dı́gitos:
>> x=22.3829
x =
22.382899999999999
>> xs = 22.3830593
xs =
22.383059299999999
>> Ex = abs(x-xs)
Ex =
1.592999999999734e-04
>> Rx = Ex / abs(x)
Rx =
7.117040240539580e-06
>>
>> y = 22.3610
y =
5
22.361000000000001
>> ys = 22.360791
ys =
22.360790999999999
>> Ey = abs(y-ys)
Ey =
2.090000000016801e-04
>> Ry = Ey/abs(y)
Ry =
9.346630293890258e-06
Vemos que los dos errores relativos son del orden de 10−6 y menores a 10−5 . Se dice que
x̂ y ŷ tienen 5 cifras significativas correctas. ¿Qué ocurre cuando restamos x̂ y ŷ? Veamos
Si le interesa leer un poco más sobre los temas de este capı́tulo, le recomendamos
la sección 1.3 del libro de John H. Mathews, Kurtis D. Fink, Métodos numéricos con
MATLAB, Pearson Educación, Prentice-Hall, 2000.
6
Capı́tulo 2
Para resolver este problema recordemos que el volumen de una esfera de radio r es
4πr3 /3, por lo que la masa de la esfera es Me = 4πρr3 /3. Si una esfera se sumerge d
centı́metros, entonces el volumen del lı́quido desalojado es
Z d
πd2 (3r − d)
Vd = π(r2 − (x − r)2 ) dx [cm3 ] = [cm3 ].
0 3
Según el principio de Arquı́medes, todo cuerpo que se sumerge en un lı́quido recibe un
empuje de abajo hacia arriba igual al peso del volumen del lı́quido desalojado. Luego, la
esfera flotará cuando se sumerja una profundidad d que haga que el peso del volumen
del lı́quido desalojado sea igual al peso de la esfera. Como la densidad del agua pura es
1 g/cm3 , el peso del volumen del lı́quido desalojado al sumergirse d centı́metros será
πd2 (3r − d)
[g]
3
y el peso de la esfera es
4πρr3 /3 [g].
Luego, la ecuación que queremos resolver es:
πd2 (3r − d)
= 4πρr3 /3
3
que equivale a
4r3 ρ − 3rd2 + d3 = 0.
Notemos que r y ρ se conocen, y que la incógnita es d. Igualmente conviene hacer todo
este razonamiento en abstracto, por si uno quiere resolverlo para otra densidad u otro
7
radio. Para ver si el problema tiene solución, es decir, si existe d entre 0 y 20(= 2r) tal
que 4r3 ρ − 3rd2 + d3 = 0 podemos graficar la función f (x) = 4r3 ρ − 3rx2 + x3 sobre
ese intervalo; aquı́ cambiamos d por x porque estamos acostumbrados a usar x para la
incógnita. En MATLAB se hace ası́:
>> r = 10
r =
10
>> rho = 0.638
rho =
0.6380
>> f = @(x) 4*r^3*rho - 3*r*x.^2 + x.^3
f =
@(x)4*r^3*rho-3*r*x.^2+x.^3
>> x=linspace(0,2*r,101);
>> plot(x,f(x),’*-’)
>> grid
>> print -dpng esferasehunde
Observación 4 (sobre MATLAB). Algunos comentarios sobre las lı́neas de MATLAB que
acabamos de ver:
La lı́nea plot(x,f(x),’*-’)
genera el gráfico que se ve a la derecha. Para ello, se arma una tabla de valores con los x
ya calculados y con los valores de f (x). Para ello, en la definición de f (x) es importante
que aparezca x.^2 y x.^3, de manera que MATLAB eleve al cuadrado y al cubo cada
componente de x. La opción ’*-’ indica que por cada punto de la tabla de valores se
debe dibujar un * y que los mismos deben ser unidos por una lı́nea -.
Finalmente la lı́nea print -dpng esferaseunde (opcional) sirve para generar el fichero
gráfico esferaseunde.png que puede ser insertado en un documento de Word, OpenOf-
fice, LaTeX, Powerpoint, etc. Nosotros la usamos para generar el gráfico que fue incluido
en este apunte.
Vemos en el gráfico que la gráfica corta el eje de las x en un punto cercano a x = 12.
Por lo tanto, dˆ = 12 cm es una primera aproximación a la solución del problema.
El objetivo de este capı́tulo es aprender sobre métodos iterativos para resolver ecua-
ciones donde no se pueda despejar la incógnita. En general, estos métodos parten de una
aproximación inicial p0 y construyen una sucesión p1 , p2 , . . . de aproximaciones cada vez
más precisas de la solución exacta p.
8
2.1. Métodos de punto fijo
El primer método que veremos es el método en que hay una fórmula, representada
por una función g que, dada una aproximación pn nos dice cómo calcular pn+1 . Es decir,
si comenzamos con una primera aproximación p0 ,
p1 = g(p0 )
p2 = g(p1 )
p3 = g(p2 )
p4 = g(p3 )
..
.
pn = g(pn−1 )
pn+1 = g(pn )
..
.
pn = 33(0.8)n → 0, cuando n → ∞.
9
Ejemplo 7. Consideremos g(x) = cos(x) y p0 = 1, hagamos este ejemplo en MATLAB:
>> p=1 .
p = .
1 .
>> p = cos(p) >> p = cos(p)
p = p =
0.5403 0.7393
>> p = cos(p) >> p = cos(p)
p = p =
0.8576 0.7389
>> p = cos(p) >> p = cos(p)
p = p =
0.6543 0.7392
>> p = cos(p) >> p = cos(p)
p = p =
0.7935 0.7390
>> p = cos(p) >> p = cos(p)
p = p =
0.7014 0.7391
. >> p = cos(p)
. p =
. 0.7391
Nos preguntamos entonces, ¿a qué converge una iteración de la forma pn+1 = g(pn )?
Observemos nuevamente lo que quiere decir iteración de la forma pn+1 = g(pn ):
p1 = g(p0 )
p2 = g(p1 )
..
.
pn = g(pn−1 )
pn+1 = g(pn )
↓ ↓
p = g(p)
Vemos que si pn → p entonces lo que está a la izquierda del signo igual converge a p.
Ahora bien, si g es una función continua y pn → p, entonces g(pn ) → g(p). Es decir, lo
que está a la derecha del signo igual, converge a g(p).
Por lo tanto, si la iteración de punto fijo pn+1 = g(pn ) genera una sucesión convergente
a p, resulta que
p = g(p).
10
Definición 8 (Punto fijo). Un pun-
to fijo de una función g es un núme-
ro p tal que p = g(p). Geométrica-
mente, p es un punto fijo si es la
abscisa (y también la ordenada) de
un punto intersección entre la gráfi-
ca de y = x y la de y = g(x).
Teorema 9. Supongamos que g es una función continua y {pn }∞ n=0 es una sucesión
generada por la iteración de punto fijo (es decir pn+1 = g(pn )). Si pn → p cuando n → ∞,
entonces p es un punto fijo de g, es decir p = g(p).
Ahora nos preguntamos: ¿Qué hace falta para poder asegurar que g tiene un punto
fijo en un intervalo [a, b]? El teorema siguiente da una respuesta a esta pregunta.
Observación 11. El Teorema 9 nos dice que el método de punto fijo con una función g
converge a un p tal que p = g(p). Por lo tanto, si queremos resolver una ecuación por el
método de punto fijo, debemos transformarla en una ecuación de la forma x = g(x). Por
ejemplo, si queremos resolver
x4 + x − sen(x) + log(1 + x2 ) = 8
x = 8 − x4 + sen(x) − log(1 + x2 ).
11
Ejemplo 12. Consideremos ahora la función g(x) = e−x en el intervalo [0, 1]. ¿Tiene un
punto fijo? Sı́, pues g(0) = 1 > 0 y g(1) = e−1 ≈ 0.3679 < 1. Veamos qué ocurre al iterar
pn+1 = e−pn comenzando con p0 = 0.1.
>> p=0.1
p =
0.1000
>> p = exp(-p)
p =
0.9048
>> p = exp(-p)
p =
0.4046
>> p = exp(-p)
p =
0.6672
>> p = exp(-p)
p =
0.5131
>> p = exp(-p)
p =
0.5986
>> p = exp(-p)
p =
0.5496
.
.
.
p =
0.5670
>> p = exp(-p)
p =
0.5673
>> p = exp(-p)
p =
0.5671
>> p = exp(-p)
p =
0.5672
>> p = exp(-p)
p =
0.5671
>> p = exp(-p)
p =
0.5672
>> p = exp(-p)
p =
0.5671
>> p = exp(-p)
p =
0.5671
12
La iteración converge a p = 0.5671 que cumple p = e−p .
Nos preguntamos ahora, ¿en qué casos podemos asegurar que la iteración de punto
fijo converge?
Analicemos los siguientes casos gráficamente:
13
14
15
En base a estos ejemplos vemos que la convergencia tiene que ver con el valor de la
derivada de g cerca del punto fijo. El resultado que vale es el siguiente:
Teorema 15 (Teorema del punto fijo). Si g es C 1 en un entorno [p − δ, p + δ] de un
punto fijo p de g y
luego
|p − p1 | = |g(p) − g(p0 )| = |g 0 (z1 )| |p − p0 | ≤ K|p − p0 | ≤ Kδ,
y por lo tanto p1 ∈ [p − δ, p + δ]. De la misma manera, p2 , p3 , · · · ∈ [p − δ, p + δ].
Además, por ese razonamiento, vemos que
|p − p1 | ≤ K|p − p0 |,
y de la misma manera
|p − p2 | ≤ K|p − p1 |
|p − p3 | ≤ K|p − p2 |
|p − p4 | ≤ K|p − p3 |
..
.
|p − pn | ≤ K|p − pn−1 |
|p − pn+1 | ≤ K|p − pn |
..
.
|p − p2 | ≤ K|p − p1 | ≤ K K|p − p0 | = K 2 |p − p0 |
|p − p3 | ≤ K|p − p2 | ≤ K K 2 |p − p0 | = K 3 |p − p0 |
|p − p4 | ≤ K|p − p3 | ≤ K K 3 |p − p0 | = K 4 |p − p0 |
16
y en general
|p − pn | ≤ K n |p − p0 |
Como 0 ≤ K < 1, resulta que K n → 0 cuando n → ∞ y entonces |p − pn | → 0.
En los libros de cálculo numérico suele aparecer esta otra versión del teorema del
punto fijo.
Teorema 17 (Punto fijo, otra versión). Sea g una función C 1 en un intervalo [a, b] que
cumple
g(x) ∈ [a, b], para todo x ∈ [a, b].
Si existe K < 1 tal que |g 0 (x)| ≤ K para todo x ∈ [a, b], entonces:
Comenzando a partir de cualquier p0 ∈ [a, b], la iteración de punto fijo pn+1 = g(pn )
converge p y más aún,
|p − pn | ≤ K n |p − p0 |. (2.1)
|p − p0 | ≤ |b − a|.
|p − pn | ≤ K n |b − a|.
p − pn = p − pn+1 + pn+1 − pn
= g(p) − g(pn ) + pn+1 − pn
= g 0 (zn ) p − pn + pn+1 − pn
p − pn − g 0 (zn ) p − pn = pn+1 − pn .
17
Es decir
1 − g 0 (zn ) p − pn = pn+1 − pn ,
que implica
1
p − pn = 0
pn+1 − pn .
1 − g (zn )
= g(p) − g(pn ) = g 0 (zn ) p − pn resulta
Como p − pn+1
g 0 (zn ) g 0 (p)
p − pn+1 = p n+1 − p n ≈ p n+1 − p n ,
1 − g 0 (zn ) 1 − g 0 (p)
0
0
g (p)
Si |g (p)| ≤ 1/2 resulta
≤ 1 y luego
1 − g 0 (p)
A continuación mostramos una implementación del método de punto fijo para resolver
la ecuación
cos x
x= .
2
Es importante notar que no se almacenan todas las iteraciones sino las últimas dos.
La última en la variable p y la anterior a la última en p0. El algoritmo se detiene entonces
cuando abs(p-p0) ≤ tol . Para asegurarnos que esto ocurra, lo que está dentro del
bloque while se ejecuta mientras abs(p-p0) > tol .
El control por el número máximo de iteraciones se realiza porque puede ocurrir que
estemos ante un caso en que el método de punto fijo no converge, y entonces hacemos
que finalice cuando se supera este número.
Al finalizar el bloque while, hay un bloque if. Este bloque pregunta: si abs(p-p0) > tol,
eso quiere decir que el while terminó porque se alcanzó el número máximo de iteraciones.
En ese caso se informa Se alcanzo el nro maximo de iteraciones. En caso contrario
(else) quiere decir que el bloque while finalizó porque se alcanzó la tolerancia desea-
da, es decir no se cumple más que abs(p-p0) > tol, luego abs(p-p0) ≤ tol , que
equivale a |pn+1 − pn | ≤ tol.
18
% iteracion inicial
p0 = 1; % elegir bien el p0 inicial <--
% % Comienza la resolucion
%p almacena la ultima iteracion
% p0 la anterior
p = g(p0);
contador = 1;
while ((abs(p-p0) > tol) && (contador < maxiter))
p0 = p;
p = g(p0);
contador = contador + 1;
end
f (x) = 0. (2.3)
Si p es un número tal que f (p) = 0 se dice que p es un cero de f o que p es una raı́z de
la ecuación (2.3).
Antes de ver el método, recordamos el teorema del valor intermedio, que es la base
que da la idea al método de bisección.
Corolario 21. Si f : [a, b] → R es una función continua y f (a) y f (b) tienen signo
opuesto (f (a) · f (b) < 0), entonces existe p ∈ [a, b] tal que f (p) = 0.
19
La idea del método de bisección es comenzar con un intervalo [a, b] tal que f (a) y f (b)
tienen signo opuesto, es decir, tenemos localizada una raı́z en el intervalo [a, b]. Luego,
ir reduciendo sistemáticamente la longitud del intervalo donde estamos seguros que hay
una raı́z, hasta que tengamos una tolerancia aceptable, es decir, un intervalo tan pequeño
donde se encuentra la raı́z que es satisfactorio para el problema en cuestión. La manera
de subdividir el intervalo consiste en tomar el punto medio del intervalo
a+b
m= ,
2
y luego analizar las tres posibilidades que pueden darse
(1) Si f (a) y f (m) tienen signos opuestos, entonces hay un cero en [a, m].
(2) Si f (m) y f (b) tienen signos opuestos, entonces hay un cero en [m, b].
···
a0 ≤ a1 ≤ · · · ≤ an ≤ · · · ≤ b n ≤ · · · ≤ b 1 ≤ b 0 .
20
• Definimos [an+1 , bn+1 ] = [an , mn ] (o an+1 = an y bn+1 = mn )
Observamos que con este procedimiento, en cada paso, la longitud del intervalo se
reduce a la mitad, es decir
|bn − an |
|bn+1 − an+1 | = .
2
Por lo tanto, después de n pasos, la longitud del intervalo es
n
|b0 − a0 | 1
|bn − an | = = |b0 − a0 |.
2n 2
Como mn está en el centro del intervalo [an , bn ] resulta que la distancia a la raı́z p es
menor o igual a la mitad de la longitud del intervalo [an , bn ], es decir
n+1
|bn − an | 1
|p − mn | ≤ = |b0 − a0 |.
2 2
x sen(x) = 1
x sen(x) − 1 = 0,
es decir, f (x) = x sen(x). Notamos que f (0) = −1 < 0 y f (2) = 2 sin(2) − 1 ≈ 0.81859 >
0. Hicimos esto en OCTAVE con las siguientes lı́neas:
21
Por lo tanto, hay una solución en el intervalo [0, 2] a la que convergerá el método de
bisección.
if (f(a)*f(b) > 0)
disp(’f(a) y f(b) tienen el mismo signo!!!!!!’)
return
end
% % Resolucion
contador = 1;
I = abs(b-a)/2;
m = (a+b)/2;
while ( (I > tol) && (contador < maxiter) )
if ( f(m) == 0 ) % el doble = es para "preguntar"
break; % sale del while de una
end
if ( f(a)*f(m) < 0 )
b = m; % a queda igual a a
else % f(b)*f(m) < 0 (no queda otra)
a = m; % b queda igual a b
end
m = (a+b)/2;
contador = contador + 1;
I = I/2;
end
if (contador == maxiter)
disp(’Se alcanzo el nro maximo de iteraciones’)
else
disp(’La solucion esta almacenada en la variable m’)
m
end
22
El código del método de bisección presentado es lo que se conoce en MATLAB como un
script. Es muy útil programar algunas resoluciones como scripts, pero conviene programar
algunas como functions (funciones). Cuando programemos un método de resolución de
algún problema, que puede ser utilizado para muchos otros problemas, lo haremos usando
funciones. Estas permiten la re-utilización sin necesidad de cambiar nombres de variables,
y las variables que se usan dentro de las funciones no afectan a las ya existentes.
A continuación mostramos cómo se modifica el código anterior para transformarlo en
una function.
if (f(a)*f(b) > 0)
disp(’f(a) y f(b) tienen el mismo signo!!!!!!’)
return
end
% % Resolucion
contador = 1;
I = abs(b-a)/2;
m = (a+b)/2;
while ( (I > tol) && (contador < maxiter) )
if ( f(m) == 0 ) % el doble = es para "preguntar"
break; % sale del while de una
end
if ( f(a)*f(m) < 0 )
b = m; % a queda igual a a
else % f(b)*f(m) < 0 (no queda otra)
a = m; % b queda igual a b
end
m = (a+b)/2;
contador = contador + 1;
I = I/2;
end
if (contador == maxiter)
disp(’Se alcanzo el nro maximo de iteraciones’)
end
23
La primera lı́nea dice:
Esto quiere decir que la función bisec recibe cinco argumentos de entrada: f, a, b,
tol, maxiter.
Para usar la función bisec podemos hacer lo siguiente:
• f = h = @(x) cos(x) - x
• a = 0
• b = 1
• tol = 1e-8
• maxiter = 100
Otra observación importante, las primeras lı́neas que tienen % constituyen el help
de la función, hasta la primera lı́nea en blanco. Ası́, ocurre lo siguiente:
24
>> help bisec
’bisec’ is a function from the file /home/pmorin/Dropbox/0-Matematica-D/codigos/b
Una diferencia importante entre las funciones y los scripts es que en las funciones
las variables son locales. Esto quiere decir que no interfieren con las variables que
tengamos definidas. Por ejemplo, observemos lo que ocurre con las siguientes lı́neas:
>> f = [1 2 3];
>> I = 80;
>> x = bisec( @(x) cos(x)-x , 0, 1, 1e-8, 100)
x = 0.739085130393505
>> f
f =
1 2 3
>> I
I = 80
>> > m
error: ’m’ undefined near line 1 column 1
Script esferasehunde.m.
format compact
r = 10;
rho = 0.638;
f = @(d) 4*r^3*rho - d.^2*3*r + d.^3;
x = linspace(0,2*r,100);
y = f(x);
plot(x,y,’*-’)
grid
25
% resolvemos
p = bisec(f,0,2*r,1e-8,1000);
f (x) = 0.
26
Para hallar la fórmula de p1 en términos de p0 observamos lo siguiente:
La recta tangente a la gráfica de y = f (x) en el punto (p0 , f (p0 )) tiene pendiente
m = f 0 (p0 ). Como la recta pasa por los puntos (p0 , f (p0 )) y (p1 , 0) resulta que
0 − f (p0 )
= m = f 0 (p0 ).
p1 − p0
Despejando p1 obtenemos
f (p0 )
p1 = p0 − .
f 0 (p0 )
Ahora podemos repetir este procedimiento, y obtener p2 a partir de p1 , la fórmula resulta
f (p1 )
p2 = p1 − .
f 0 (p1 )
En general,
f (pn )
pn+1 = pn − .
f 0 (pn )
Observación 23. El método de Newton es un método de punto fijo, con función de
iteración dada por
f (x)
g(x) = x − 0 .
f (x)
Notar que x = g(x) si y sólo si f (x) = 0 (siempre que f 0 (x) 6= 0). En otras palabras, x es
solución de la ecuación original f (x) = 0 si y sólo si x es un punto fijo de g.
Para comprender cómo funciona el método de Newton, analicemos g 0 (x):
Luego, si g 0 es una función continua, resulta que existe δ > 0 tal que |g 0 (x)| < 1/2 en
[p − δ, p + δ].
Para que g 0 sea continua es suficiente que f sea C 2 y que f 0 no se anule en una
vecindad de x = p.
Recordando el Teorema 15 hemos deducido que
Teorema 24. [Convergencia del método de Newton] Sea p un cero de la función f , es
decir
f (p) = 0.
27
Si f es C 2 en una vecindad de p, f 0 (p) 6= 0 y p0 está suficientemente cerca de p, entonces
la iteración de Newton
f (pn )
pn+1 = pn − , n = 0, 1, 2, . . .
f 0 (pn )
converge a p.
Observación 25. El método de Newton es un método de punto fijo que siempre funciona
(si comenzamos suficientemente cerca). Nos da una manera de construir o elegir la función
de iteración g para que siempre resulte convergente.
Veamos un ejemplo, volvamos a resolver la ecuación
x = cos(x).
Resolvámosla primero con el método de punto fijo pn+1 = cos(pn ). Para ir contando
iteraciones lo hacemos de esta manera:
>> p = 1, cont = 0, .
p = 1 .
cont = 0 .
>> p = cos(p), cont = cont + 1 >> p = cos(p), cont = cont + 1
p = 0.54030 p = 0.73908
cont = 1 cont = 31
>> p = cos(p), cont = cont + 1 >> p = cos(p), cont = cont + 1
p = 0.85755 p = 0.73909
cont = 2 cont = 32
>> p = cos(p), cont = cont + 1 >> p = cos(p), cont = cont + 1
p = 0.65429 p = 0.73908
cont = 3 cont = 33
>> p = cos(p), cont = cont + 1 >> p = cos(p), cont = cont + 1
p = 0.79348 p = 0.73909
cont = 4 cont = 34
>> p = cos(p), cont = cont + 1 >> p = cos(p), cont = cont + 1
p = 0.70137 p = 0.73908
cont = 5 cont = 35
>> p = cos(p), cont = cont + 1 >> p = cos(p), cont = cont + 1
p = 0.76396 p = 0.73909
cont = 6 cont = 36
. >> p = cos(p), cont = cont + 1
. p = 0.73909
. cont = 37
El método de punto fijo necesitó 37 iteraciones para que la diferencia entre dos ite-
randos consecutivos fuera menor a 10−5 .
Veamos ahora qué ocurre si usamos el método de Newton. Escribimos primero la
ecuación de la forma
x − cos(x) = 0.
Es decir, f (x) = x − cos(x) y f 0 (x) = 1 + sen(x). Luego, la iteración de Newton toma la
forma:
pn − cos(pn )
pn+1 = pn − .
1 + sen(pn )
28
En OCTAVE
>> p = 1, cont = 0
p = 1
cont = 0
>> p = p - (p - cos(p))/(1 + sin(p)), cont = cont + 1
p = 0.75036
cont = 1
>> p = p - (p - cos(p))/(1 + sin(p)), cont = cont + 1
p = 0.73911
cont = 2
>> p = p - (p - cos(p))/(1 + sin(p)), cont = cont + 1
p = 0.73909
cont = 3
>> p = p - (p - cos(p))/(1 + sin(p)), cont = cont + 1
p = 0.73909
cont = 4
Newton no deja de sorprendernos. Bastaron cinco iteraciones para que el error entre
dos consecutivos sea menor que 10−16 . Impresionante. . .
Veamos lo que ocurre con el error
29
n p |p − pn |
0 1.00000000000000000000 0.26091486678483932771
1 0.75036386784024389218 0.01127873462508321989
2 0.73911289091136167517 0.00002775769620100288
3 0.73908513338528392111 0.00000000017012324882
4 0.73908513321516056127 0.00000000000000011102
5 0.73908513321516067229 0.00000000000000000000
Vemos que la cantidad de ceros después del punto decimal en el error es, en cada
iteración, aproximadamente el doble que en la iteración anterior.
Los números de esa tabla fueron generados con estas lı́neas de código OCTAVE
pp = 0.73908513321516067229;
p = 1;
for i=0:5
fprintf(’%d %22.20f %22.20f \n’,i, p, abs(p-pp));
p = p - (p - cos(p))/(1 + sin(p));
end
Para entender por qué es tan rápido el método de Newton, definiremos a continuación
una manera de estudiar el orden de convergencia de los métodos iterativos.
30
Vemos que en el caso A = 0.1 el error es en cada paso un diez porciento √ del error en el
caso anterior. Se dice que se gana un dı́gito por iteración. En el caso A = 0.1 ≈ 0.316 se
gana un dı́gito cada dos iteraciones. Lo importante es que en los métodos de convergencia
lineal se gana un dı́gito cada un cierto número de iteraciones, pero esto se mantiene fijo
a lo largo de toda la iteración.
Veamos qué ocurre con un caso de convergencia cuadrática, es decir, con R = 2. Para
simplificar la presentación consideramos A = 1.
A = 1, R = 2, |p − pn+1 | ≈ |p − pn |2
n |p − pn |
0 10−1 = 0.1
1 10−2 = 0.01
2 10−4 = 0.0001
3 10−8 = 0.00000001
4 10−16 = 0.0000000000000001
5 10−32 = 0.00000000000000000000000000000001
6 10−64 = . . .
|f 00 (p)|
|p − pn+1 | ≈ |p − pn |2 .
2|f 0 (p)|
31
y por el teorema de Taylor
(x − pn )2 00
f (x) = f (pn ) + (x − pn )f 0 (pn ) + f (z)
2
para algún z entre x y pn . Esta fórmula nos dice que la ecuación de la recta tangente nos
da una aproximación de la función f (x) con un error
(x − pn )2 00
f (z).‘
2
2
Este error es menor a M2 (x−p
2
n)
si M2 es una cota para el valor absoluto de f 00 . Como ese
error es cuadrático en términos de la distancia a pn , el método de Newton resulta con
convergencia cuadrática.
A continuación presentamos el código del método de Newton.
Función newton.m.
function p = newton(f,fprima,p0,tol,maxiter)
% p = newton(f,fprima,p0,tol,maxiter)
% metodo de Newton para ecuaciones no lineales
% f: funcion escalar. Ecuacion a resolver
% f(x) = 0
% fprima: derivada de f
% p0: aproximacion inicial
% tol: tolerancia para el error absoluto
% maxiter: maximo numero de iteraciones permitido
contador=1;
p=p0-f(p0)/fprima(p0);
if(contador>maxiter)
disp(’se supero el maximo de iteraciones’)
end
32
La idea es comenzar entonces con p0 y p1 cercanos a p, y definir p2 como la abscisa del
punto intersección entre el eje x y la recta que pasa por (p0 , f (p0 )) y por (p1 , f (p1 )). Para
obtener la fórmula de p2 en términos de p0 y p1 observamos primero que la pendiente de
la recta secante es
f (p1 ) − f (p0 )
m= .
p 1 − p0
Luego, usando que (p1 , f (p1 )) y (p2 , 0) pertenecen a esa recta, resulta
f (p1 ) − 0
= m.
p1 − p2
En consecuencia,
f (p1 ) − 0 f (p1 ) − f (p0 )
=
p1 − p2 p1 − p0
y despejando p2 obtenemos
p 1 − p0
p2 = p1 − f (p1 ) .
f (p1 ) − f (p0 )
33
Teorema 29. [Convergencia del método de la secante] Sea p un cero de la función f , es
decir
f (p) = 0.
pn − pn−1
pn+1 = pn − f (pn ) , n = 1, 2, . . .
f (pn ) − f (pn−1 )
√
1+ 5
converge a p. El orden del método es R = 2
≈ 1.61803 . . . (número de oro o razón
áurea).
Función secante.m.
function p = secante(f,p0,p1,tol,maxiter)
% p = secante(f,p0,p1,tol,maxiter)
% metodo de la secante para ecuaciones no lineales
% f: funcion escalar. Ecuacion a resolver
% f(x) = 0
% p0, p1: aproximaciones iniciales
% tol: tolerancia para el error absoluto
% maxiter: maximo numero de iteraciones permitido
contador=1;
f0 = f(p0);
f1 = f(p1);
p = p1 - f1 * (p1-p0) / (f1-f0);
if(contador>maxiter)
disp(’se supero el maximo de iteraciones’)
end
34
2.6. Métodos iterativos para sistemas de ecuaciones
Consideremos el siguiente sistema de dos ecuaciones no lineales con dos incógnitas:
x2 − 2x − y − 0.5 = 0,
x2 + 4y 2 − 4 = 0.
close all
contour(x,y,z1,[0 0],’b’)
hold on
contour(x,y,z2,[0 0],’r’)
hold off
grid
Con ayuda del comando help investigue qué hacen las instrucciones meshgrid y
contour para entender qué está graficando y cómo harı́a para graficar curvas en otro
problema.
Vemos a partir del gráfico que hay dos soluciones. Se ve que podrı́a haber 3, 4, o
ninguna, dependiendo de la posición relativa de la parábola con respecto a la elipse.
Para otros sistemas de ecuaciones puede ocurrir que haya un número n a priori des-
conocido de soluciones, o también infinitas.
Recordemos que para sistemas lineales hay sólo tres opciones:
No hay solución.
35
punto fijo. Para ello, podemos despejar un x de la primera ecuación y un y de la segunda,
por ejemplo:
x2 − y + 0.5 −x2 − 4y 2 + 8y + 4
x= , y= .
2 8
Ahora reescribimos estas dos ecuaciones como una ecuación vectorial
" x2 −y+0.5 #
x
= −x2 −4y22 +8y+4 .
y 8
0
Comenzamos ahora las iteraciones de punto fijo a partir de p = que está cerca de
1
uno de las soluciones del sistema de ecuaciones. Lo hacemos en OCTAVE:
>> g = @(p) [ (p(1)^2-p(2)+0.5)/2 >> p = g(p)
(-p(1)^2-4*p(2)^2+8*p(2)+4)/8 ]; p =
>> p = [0;1] -0.22219
p = 0.99380
0 >> p = g(p)
1 p =
>> p = g(p) -0.22222
p = 0.99381
-0.25000 >> p = g(p)
1.00000 p =
>> p = g(p) -0.22221
p = 0.99381
-0.21875 >> p = g(p)
0.99219 p =
>> p = g(p) -0.22221
p = 0.99381
-0.22217 >> p = g(p)
0.99399 p =
>> p = g(p) -0.22221
p = 0.99381
-0.22231
0.99381
2
Veamos ahora qué ocurre cuando comenzamos con que está cerca de la otra
0
solución.
36
>> p = [2;0] >> p = g(p)
p = p =
2 1011.99
0 -392.60
>> p = g(p) >> p = g(p)
p = p =
2.25000 5.1226e+05
0.00000 -2.0548e+05
>> p = g(p) >> p = g(p)
p = p =
2.78125 1.3121e+11
-0.13281 -5.3912e+10
>> p = g(p) >> p = g(p)
p = p =
4.18408 8.6076e+21
-0.60855 -3.6052e+21
>> p = g(p) >> p = g(p)
p = p =
9.3075 3.7046e+43
-2.4820 -1.5760e+43
>> p = g(p) >> p = g(p)
p = p =
44.806 6.8619e+86
-15.891 -2.9574e+86
Vemos que en este caso la iteración no
converge.
Veamos qué ocurre si comenzamos
1.9
más cerca, por ejemplo, a partir de p = :
0.31
octave:54> p=[1.9;0.31] .
p = .
1.90000 .
0.31000 octave:71> p = g(p)
octave:55> p = g(p) p =
p = -0.22192
1.90000 0.99378
0.31070 octave:72> p = g(p)
octave:56> p = g(p) p =
p = -0.22227
1.89965 0.99382
0.31118 octave:73> p = g(p)
octave:57> p = g(p) p =
p = -0.22221
1.89874 0.99381
0.31168 octave:74> p = g(p)
octave:58> p = g(p) p =
. -0.22221
. 0.99381
.
En este caso, la iteración converge a la primera solución que habı́amos encontrado.
Claramente, el punto que queremos hallar, no es un punto de atracción de esa iteración
de punto fijo.
37
En el caso de una ecuación con una incógnita, hemos visto que la iteración de punto
fijo pn+1 = g(pn ) converge localmente al punto fijo p si |g 0 (p)| < 1. Se dice que un método
converge localmente si lo hace cuando comenzamos suficientemente cerca.
En el caso de sistemas de ecuaciones, el rol la derivada g 0 lo juega la matriz Jacobiana:
Si
2 2 x1 g1 (x1 , x2 )
G:R →R , G = ,
x2 g2 (x1 , x2 )
su matriz Jacobiana o su diferencial es la matriz
∂g ∂g1
1
(x1 , x2 ) (x1 , x2 )
0 x 1 ∂x 1 ∂x2
G = ∂g
x2 2 ∂g2
(x1 , x2 ) (x1 , x2 )
∂x1 ∂x2
Para recordar, la primera fila de G0 tiene todas las derivadas de g1 y la segunda fila de
G0 tiene todas las derivadas de g2 .
En general, si G : RN → RN y
g1 (x) x1
g2 (x1 ) x2
G(x) = .. , con x = .. ,
. .
gN (x1 ) xN
p = G(p).
Si G es C 1 en un entorno de p y
N
0 0
X ∂gi
kG (p)k1 < 1, con kG (p)k1 = máx ∂x (p) ,
1≤i≤N j
j=1
pn+1 = G(pn )
converge localmente. Es decir, hay un entorno alrededor de p tal que si p0 está en ese
entorno, resulta que pn → p.
38
Observación 31. En general es difı́cil de verificar que kG0 (p)k1 < 1, sobre todo porque
no se conoce p. El teorema recién enunciado es de carácter teórico.
Ahora pensamos, si estamos dispuestos a calcular todas esas derivadas, ¿por qué no
usamos un método de Newton?
f1 (x1 , x2 , . . . , xN ) = 0
f2 (x1 , x2 , . . . , xN ) = 0
..
.
fN (x1 , x2 , . . . , xN ) = 0,
F (x) = 0.
39
F (p0 ) es un vector columna de N componentes, o una matriz de N × 1.
p1 = p0 + δp.
F 0 (pn ) δp = − F (p0 ) .
| {z } |{z} | {z }
conocido incógnita conocido
Función newtonvec.m.
function p= newtonvec(F,Fprima,p0,tol,maxiter)
% p = newtonvec(F,Fprima,p0,tol,maxiter)
% metodo de Newton para sistemas
% F recibe un vector columna y devuelve un vector columna
% Fprima: recibe un vector columna y devuelve una matriz,
% la matriz Jacobiana de F
% p0: aproximacion unicial (vector columna)
% tol: tolerancia para el error absoluto
% maxiter: maximo numero de iteraciones permitido
contador = 1;
deltap = - Fprima(p0) \ F(p0);
p = p0 + deltap;
40
Veamos cómo se utiliza esta función para hallar la solución del siguiente sistema de
ecuaciones:
x2 y + xy 3 = 9
3x2 y − y 3 = 4.
Primero las escribimos con lado derecho igual a cero:
x2 y + xy 3 − 9 = 0
3x2 y − y 3 − 4 = 0.
Luego, cambiamos x por x1 y y por y2 :
x21 x2 + x1 x32 − 9 = 0
3x21 x2 − x32 − 4 = 0.
Luego, definimos F : R2 → R2 ,
x21 x2 + x1 x32 − 9
F (x) =
3x21 x2 − x32 − 4
0
y el sistema es equivalente a resolver F (x) = . El Jacobiano de F es
0
0 2x1 x2 + x32 x21 + 3x1 x22
F (x) =
6x1 x2 3x21 − 3x22
1
Comenzaremos con la aproximación inicial p0 = . Ya estamos en condiciones de
1
hacerlo en OCTAVE. Como hay que hacer un par de definiciones medio largas, conviene
poner todo en un script en lugar de tipear en la lı́nea de comandos.
Script pruebanewtonvec.m.
% pruebanewtonvec.m
% prueba Newton vectorial
% Resolvemos el problema
% x1^2 + x1 x2^3 - 9 = 0
% 3 x1^2 x2 - x2^3 - 4 = 0
% luego su Jacobiano
% en el Jacobiano, cada fila tiene las derivadas de cada
% fila de F, con respecto a cada variable.
DF = @(x) [ 2*x(1)+x(2)^3 3*x(1)*x(2)^2
6*x(1)*x(2) 3*x(1)^2-3*x(2)^2]
p = newtonvec(F,DF,[1;1],1e-12,100)
F(p)
41
Como el archivo se llama pruebanewtonvec.m para ejecutarlo escribimos en la lı́nea
de comandos lo siguiente:
>> pruebanewtonvec
F =
@(x) [x(1) ^ 2 + x(1) * x(2) ^ 3 - 9; 3 * x(1) ^ 2 * x(2) - x(2) ^ 3 - 4]
DF =
@(x) [2 * x(1) + x(2) ^ 3, 3 * x(1) * x(2) ^ 2; 6 * x(1) * x(2), 3 * x(1) ^ 2 - 3 * x(
p =
1.3364
1.7542
ans =
0
0
>> format long
>> p
p =
1.33635537721717
1.75423519765170
2.7. Ejercicios
2.1. (a) Graficar la función g(x) = 4.8 − |x − 3.5|2.5 sobre el intervalo [1, 5].
(2 + sen(2x))2
2.2. (a) Graficar la función g(x) = 2
y su derivada g 0 (x) en gráficos sepa-
8(1 + cos (x))
rados. Observar que hay un único punto fijo.
(b) Determinar, a partir del gráfico, un valor de K tal que |g 0 (x)| ≤ K para todo x ∈
[0, 1].
42
2.3. Consideremos la ecuación
x4 + 2x2 − x − 3 = 0.
(a) Verificar que α ∈ [1, 32 ] es solución de esa ecuación si y solo si α = gi (α) para cada
una de las funciones gi (x) siguientes. Es decir, α ∈ [1, 23 ] es solución de esa ecuación
si y solo si α es un punto fijo de gi .
1/2
x + 3 − x4
2 1/4
g1 (x) = (3 + x − 2x ) g2 (x) =
2
1/2
3x4 + 2x2 + 3
x+3
g3 (x) = g4 (x) =
x2 + 2 4x3 + 4x − 1
(b) Verificar con cuál de ellas se obtiene una convergencia más rápida a la solución
comenzando de p0 = 1.
2.4. Baje el código de la función p = newton(f, fprima, x0, tol, maxiter) del apun-
te y guárdelo en su carpeta de trabajo. Verifique que comprende su uso y cómo fue
programado.
Verifique el funcionamiento probando con un problema del cual se conozca la solución.
Por ejemplo, resolver x2 = 25. Verifique que el error decrece cuadráticamente.
3
2.5. El polinomio p(x) = −x3 + 200 x2 + x + 14 tiene tres ceros entre −2 y 2. Utilice
el método de Newton para hallarlas con un error absoluto menor a 10−12 . Para ello,
grafique el polinomio y determine a partir del gráfico tres diferentes valores iniciales p0
que aseguren la convergencia a cada cero del polinomio.
2.7. Baje el código de la función p = secante(f, p0, p1, tol, maxiter) del apunte
y guárdelo en su carpeta de trabajo. Verifique que comprende su uso y cómo fue progra-
mado.
Verifique el funcionamiento probando con un problema del cual se conozca la solución.
Por ejemplo, resolver x3 = 27.
2.8. La ecuación
ln(1 + x4 )
x
exp 2 =
x + cos(x) 2
tiene una solución en el intervalo [0, 5]. Calcúlela con un error absoluto menor a 10−12
utilizando el método de la secante. Determine p0 y p1 adecuados a partir de un gráfico.
43
2.9. Se sabe que la ecuación !
3
ln(3 + x )
x = ln 2 sen(x2 )
2+ 1+x2
tiene una solución en [0, 1]. Hállela con un error relativo menor que 10−12 .
2.10. Para cada uno de los siguientes sistemas, determinar cuántas soluciones tiene en el
rectángulo indicado, a partir de un gráfico:
(a) x2 + xy 3 = 9, 3x2 y − y 3 = 4 en el cuadrado [−4, 4] × [−4, 4].
(b) 3w2 − z 2 = 0, 3wz 2 − w3 = 1, en el rectángulo [0, 3] × [0, 2].
2.11. Escribir una función p = newtonvec(F, JF, p0, tol, maxiter) que reciba co-
mo argumentos:
f: una función vectorial que reciba un vector x y devuelva un vector y = f (x)
(ambos vectores columna); se quiere resolver F (x) = 0;
JF: la función Jacobiano de F; dado x devuelve una matriz que en la posición (i, j)
∂fi
contiene ∂x j
2.14. Resolver el siguiente sistema de ecuaciones con una precisión para el error absoluto
de 10−6 , utilizando el método de Newton:
p
6x − 2 cos(yz) = 1, 9y + x2 + sen(z) + 1.06 + 0.9 = 0, 60z + 3e−xy = 3 − 10π.
Comenzar las iteraciones con x = 0, y = 0, z = 0.
44
2.15. En un sistema cerrado tienen lugar las siguientes reacciones quı́micas reversibles:
2A + B C (1),
A+D C (2).
Si le interesa leer un poco más sobre los temas de este capı́tulo, le recomendamos las
secciones 2.1 a 2.4 y la sección 3.7 del libro de John H. Mathews, Kurtis D. Fink, Métodos
numéricos con MATLAB, Pearson Educación, Prentice-Hall, 2000.
45
Capı́tulo 3
Aproximación e Interpolación
Polinomial
En este caso, lo que se hace, es aproximar la función por funciones sencillas que se sepan
integrar, y luego se integran esas funciones sencillas.
Para los matemáticos, las funciones más sencillas son los polinomios. Algunas razones
por las que se consideran sencillas son las siguientes:
(1) Son fáciles de representar, pues basta un número finito de coeficientes. En OCTAVE,
un vector p de n componentes representa el siguiente polinomio de grado n − 1:
>> p = [6 -3 0 1];
>> x = linspace(-3,3,300);
>> y = polyval(p,x);
>> plot(x,y)
46
Resumiendo, para OCTAVE, p es un vector de 4 componentes. Cuando lo evaluamos
con polyval, lo interpreta como un polinomio de grado 3. Pero para que lo interprete
de esa manera, hay que evaluarlo con polyval. Se ruega al lector entender lo que
ocurre al hacer las siguientes lı́neas:
>> p(3)
ans = 0
>> polyval(p,3)
ans = 136
>> p(0)
error: subscript indices must be either positive integers less than 2^31 or logica
>> polyval(p,0)
ans = 1
>> p(4)
ans = 1
>> polyval(p,4)
ans = 337
>> p(0.25)
error: subscript indices must be either positive integers less than 2^31 or logica
>> polyval(p,0.25)
ans = 0.90625
d
p1 xn−1 + p2 xn−2 + · · · + pn−1 x + pn
dx
= (n − 1)p1 xn−2 + (n − 2)p2 xn−3 + · · · + 3pn−3 x2 + 2pn−2 x + pn−1 .
y = ax + b
47
cuya gráfica pasa por los puntos.
Si tengo tres puntos con diferente abscisa, por ejemplo
y = ax2 + bx + c.
Veamos cómo se calcula dicha parábola. Para que la misma pase por los puntos debe
cumplir
Nos queda entonces un sistema de 3 ecuaciones con 3 incógnitas (las incógnitas son a, b,
c):
12 a + 1b + c = 8
(−1)2 a + (−1)b + c = 9
22 a + 2b + c = 15.
>> X = [1 -1 2]’;
>> Y = [8 9 15]’;
>> matriz = [X.^2 X ones(size(X))]
matriz =
1 1 1
1 -1 1
4 2 1
>> p = matriz\Y
p =
2.50000
-0.50000
6.00000
48
Ası́, a = 2.5, b = −0.5 y c = 6. Para graficar el polinomio en el intervalo [−2, 3] y los
puntos dados hacemos lo siguiente:
>> xx = linspace(-2,3,200);
>> yy = polyval(p,xx);
>> plot(xx,yy,X,Y,’*’)
Ahora nos preguntamos qué ocurre cuando tenemos n puntos. La respuesta está en el
siguiente teorema.
con abscisas todas distintas, existe un único polinomio p de grado n − 1 cuya gráfica pasa
por los puntos, es decir
p(x1 ) = y1
p(x2 ) = y2
.. ⇐⇒ p(xi ) = yi , i = 1, 2, . . . , n.
.
p(xn ) = yn
Para demostrar este polinomio recordamos que el Teorema Fundamental del Álgebra
dice:
Sin entrar en detalles sobre lo que son los números complejos, deducimos de ese
enunciado que
p(xi ) = 0, i = 1, 2, . . . , n,
p1 = p2 = · · · = pn−1 = pn = 0.
Demostración del Teorema 33. Queremos hallar el polinomio, ası́ que proponemos
49
El mismo debe cumplir las siguiente n ecuaciones
p(x1 ) = y1
p(x2 ) = y2
..
.
p(xn ) = yn .
Es decir
p1 xn−1
1 + p2 xn−2
1 + · · · + pn−1 x1 + pn = y1
p1 xn−1
2 + p2 xn−2
2 + · · · + pn−1 x2 + pn = y2
..
.
n−1 n−2
p1 xn + p2 xn + · · · + pn−1 xn + pn = yn .
Recordando que las abscisas x1 , x2 , . . . , xn son conocidas y las ordenadas y1 , y2 , . . . , yn
también, las incógnitas son los coeficientes p1 , p2 , . . . , pn del polinomio. Para encontrarlos
debemos resolver ese sistema de ecuaciones que en forma matricial se escribe ası́:
n−1 n−2 p1 y1
x1 x1 . . . x1 1
xn−1 xn−2 . . . x2 1 p2 y2
2 2 .. ..
.. .. .. .. . = . . (3.1)
. . . .
pn−1 yn−1
xn−1
n xn−2
n . . . xn 1
pn yn
Observamos nuevamente que el lado derecho del sistema de ecuaciones es el vector de
ordenadas de los puntos (coordenadas y) y las columnas de la matriz son potencias
de las abscisas de los puntos (coordenadas x). Esta matriz se conoce como matriz de
Vandermonde correspondiente a x1 , x2 , . . . , xn .
Este sistema de ecuaciones tiene solución única si la matriz es invertible, lo cual puede
verificarse demostrando que el problema homogéneo tiene como única solución a la trivial.
Plantear el sistema homogéneo
n−1 n−2 p1
x1 x1 . . . x1 1 0
xn−1 xn−2 . . . x2 1 p2 0
2 2 .
.. .. .. .. .. = ..
. . . . .
n−1 n−2
pn−1
xn xn . . . xn 1 0
pn
equivale a buscar el polinomio que satisface p(xi ) = 0, i = 1, 2, . . . , n. Por lo que hemos
dicho antes de comenzar la demostración, esto implica que p1 = p2 = · · · = pn−1 = pn = 0,
es decir, la única solución es la trivial.
Por lo tanto, el sistema (3.1) tiene solución única, o lo que es lo mismo, existe un
único polinomio cuya gráfica pasa por los n puntos dados.
En OCTAVE, para hallar el polinomio de grado n−1 que pasa por n puntos, se utiliza
la función polyfit. Por ejemplo, para hallar el polinomio de grado 4 cuya gráfica pasa
por
(1, 2) (2, −1) (4.5, 3) (5.5, 0) (10, 3.3)
50
hacemos ası́:
>> xx = linspace(0,11,200);
>> yy = polyval(p,xx);
>> plot(xx,yy,X,Y,’*’)
Esas lı́neas también incluyen instrucciones para graficar el polinomio en el intervalo [0, 11].
51
la aproximación por cuadrados mı́nimos, que consiste en lo siguiente. Queremos hallar
el polinomio p(x) que minimice la suma de la distancia entre p(xi ) y yi al cuadrado, es
decir, que minimice la siguiente expresión:
n
2 2 2 X 2
p(x1 ) − y1 + p(x2 ) − y2 + · · · + p(xn ) − yn = p(xi ) − yi .
i=1
Observamos que
p1 x21 + p2 x1 + p3 x21 x1 1 p1
p1 x2 + p2 x2 + p3 x2 x2 1 p2
2 2
= ..
.. ..
. . .
2 2
p 1 xn + p 2 xn + p 3 xn xn 1 p n
| {z } | {z }
M p
n
!1/2
X
(M p)i − yi ]2
kM p − yk =
i=1
Para resolver este problema, recordemos que M p es una combinación lineal de las
columnas de M :
M p = p1 col1 (M ) + p2 col2 (M ) + p3 col3 (M ),
y por lo tanto, el problema consiste en encontrar el elemento del espacio columna de M
que está lo más cerca posible de y.
Recordando un poco de geometrı́a, vemos que M p es el elemento más cercano a dicho
espacio si y sólo si el error y − M p es perpendicular al espacio generado por las columnas
de M . Es decir, si y − M p es perpendicular a cada columna de M . Más precisamente, si
se cumple que
col1 (M ) · (y − M p) = 0
col2 (M ) · (y − M p) = 0
col3 (M ) · (y − M p) = 0.
52
Recordando que si u, v son vectores columna, entonces u · v = uT v resulta que M p debe
cumplir
col1 (M )T (y − M p) = 0
col2 (M )T (y − M p) = 0
col3 (M )T (y − M p) = 0.
M T M p = M T y.
Es decir, el vector p (de tres componentes) debe ser la solución al sistema de ecuaciones
MT M p = MT y
En OCTAVE, si X, Y son vectores que tienen las abscisas y las ordenadas de los puntos,
el cálculo del polinomio p se harı́a de la siguiente manera:
>> A = M’*M
(3) Hallamos los coeficientes del polinomio que aproxima en el sentido de mı́nimos cua-
drados:
>> p = A \ Y
(4) Ahora graficamos el polinomio y los puntos en un mismo gráfico. Supongamos para
ello que estamos en el intervalo [0, 6]:
>> xx = linspace(0,6,200);
>> yy = polyval(p,xx);
>> plot(X,Y,’*’,xx,yy)
Los pasos (1) a (3) de lo recién explicado pueden resumirse en OCTAVE usando
polyfit:
>> p = polyfit(X,Y,2)
53
Observación 34. Recalcamos que igualmente es importante entender cómo se resuelve
un sistema en el sentido de mı́nimos cuadrados, pues hay ocasiones en que se quieren
aproximar puntos por funciones de otra forma, por ejemplo: p(x) = a3x + b2x + c. Allı́
hay que reproducir los pasos, utilizando como matriz M la siguiente:
Observación 35. El estudio del error cometido por el polinomio hallado y la función
real que modela el comportamiento de los puntos, que han sido medidos con un cierto
error aleatorio, escapa al alcance de este curso, y corresponde al área Estadı́stica.
Función pruebapolyfit1.m.
function p = pruebapolyfit1(f,a,b,n)
% function p = pruebapolyfit1(f,a,b,n)
%
% Se interpola la funcion f en n
% puntos equiespaciados en el intervalo [a,b].
% En la figura 1 se grafica la funcion f y el polinomio.
% En la figura 2 se grafica el error.
% Ambos graficos se realizan sobre el intervalo [a,b].
%
% La funcion retorna el polinomio p en formato MATLAB.
54
% Ahora si, graficamos
figure(1)
plot(xx,fxx,’LineWidth’,2,xx,pxx,’LineWidth’,2,X,Y,’*’)
legend(’funcion original’,’polinomio interpolante’)
set(gca,’FontSize’,18)
figure(2)
plot(xx,error,’LineWidth’,2)
grid, title(’Error’)
set(gca,’FontSize’,18)
Ejemplo 36. Veamos qué ocurre cuando interpolamos la función f (x) = sen(x) sobre el
intervalo [0, 7]. Para entender mejor se recomienda que baje la función pruebapolyfit1 y
realice estas pruebas usted mismo.
Vemos en este ejemplo que el error tiende a cero a medida que aumentamos el número n
de puntos.
2
Ejemplo 37. Consideremos ahora la función f (x) = e−x sobre el intervalo [−1, 3].
55
maximo_error = 0.20698
>> pruebapolyfit1(f,-1,3,5);
maximo_error = 0.23692
>> pruebapolyfit1(f,-1,3,6);
maximo_error = 0.084354
>> pruebapolyfit1(f,-1,3,7);
maximo_error = 0.065621
>> pruebapolyfit1(f,-1,3,8);
maximo_error = 0.052012
>> pruebapolyfit1(f,-1,3,9);
maximo_error = 0.023666
>> pruebapolyfit1(f,-1,3,10);
maximo_error = 0.019472
Aquı́ también parece ser que el error tiende a cero, aunque da la impresión de que esto
ocurre más lentamente que en el ejemplo anterior.
>> pruebapolyfit1(f,-5,5,11);
maximo_error = 1.9156
>> pruebapolyfit1(f,-5,5,12);
maximo_error = 0.55418
>> pruebapolyfit1(f,-5,5,13);
56
maximo_error = 3.6554
>> pruebapolyfit1(f,-5,5,14);
maximo_error = 1.0639
>> pruebapolyfit1(f,-5,5,15);
maximo_error = 7.1908
>> pruebapolyfit1(f,-5,5,16);
maximo_error = 2.0979
>> pruebapolyfit1(f,-5,5,17);
maximo_error = 14.329
Aquı́ se ve que la cosa no anda tan bien. El error no parece tender a cero. . .
Teorema 39. Sea f una función definida en el intervalo [a, b]. Si x1 , . . . , xn son todos
puntos distintos en el intervalo [a, b], entonces existe un único polinomio p(x) de grado
menor o igual a n − 1 cuya gráfica pasa por los puntos (xi , f (xi )) para i = 0, 1, . . . , n. Es
decir, existe un único polinomio p(x) = p1 xn−1 + · · · + pn−1 x + pn tal que
p(xi ) = f (xi ), i = 1, 2, . . . , n.
Mn
|f (x) − p(x)| ≤ |x − x1 | |x − x2 | . . . |x − xn |, donde Mn = máx |f n (ξ)|.
n! ξ∈[a,b]
Corolario 40. Si no se sabe dónde están los puntos xi , entonces, para cada x ∈ [a, b]
sólo sabemos que |x − xi | ≤ |b − a|. Del teorema anterior concluimos que:
Mn
máx |f (x) − p(x)| ≤ |b − a|n .
x∈[a,b] n!
b−a
xi = a + (i − 1)∆x, con ∆x = ,
n−1
|b − a|n
|x−x1 | |x−x2 | . . . |x−xn | ≤ ∆x ∆x 2∆x 3∆x . . . (n−1)∆x = (n−1)!∆xn = (n−1)! .
(n − 1)n
Mn |b − a|n
máx |f (x) − p(x)| ≤ .
x∈[a,b] n (n − 1)n
| {z }
cota del error
57
Vale la pena observar que en ambos casos, la cota del error tiene un factor Mn que
depende de la derivada de orden n de la función que queremos aproximar, y otro factor,
|b − a|n |b − a|n |b − a|n |b − a|n
en un caso y en el otro. Tanto como tienden a cero
n! n(n − 1)n n! n(n − 1)n
cuando n tiende a infinito.
En el ejemplo 36, Mn = 1 para todo n, ası́ que
|b − a|n |b − a|n
Mn = 1 → 0, cuando n → ∞.
n(n − 1)n n(n − 1)n
2
En el ejemplo 37, resulta más pesado calcular derivadas de todos los órdenes de e−x .
Cabe aclarar que Mn → ∞ cuando n → ∞, pero aún ası́,
|b − a|n
Mn → 0, cuando n → ∞.
n(n − 1)n
Sin embargo, el valor de Mn puede tender a infinito tan rápido que los términos
completos
|b − a|n |b − a|n
Mn y Mn
n! n(n − 1)n
tiendan a infinito. Este es el caso del ejemplo 38. En ese caso, no podemos asegurar que
la aproximación sea cada vez mejor a medida que tomemos más puntos de interpolación
y polinomios de grado mayor.
Para asegurar convergencia de la aproximación, se utiliza la interpolación polinomial
a trozos.
|b − a|n
|f (x) − p(x)| ≤ Mn ,
n(n − 1)n
con Mn = máxξ∈[a,b] |f n (ξ)|. También vimos que Mn puede ser muy grande para n gran-
de. En este caso, no sirve de mucho aumentar excesivamente el número de puntos de
interpolación.
¿Qué se puede hacer para asegurar que el error tiende a cero?
Una cosa que se puede hacer es utilizar interpolación polinomial a trozos. La idea es
fijar el grado polinomial n − 1 e interpolar con ese grado (fijo) sobre L sub-intervalos.
Más precisamente. Consideremos un intervalo [a, b] y un número entero positivo L. Sub-
|b − a|
dividimos el intervalo [a, b] en L sub-intervalos de longitud h = . Llamamos y1 , y2 ,
L
. . . , yL+1 a los extremos de esos sub-intervalos, es decir
58
En el gráfico a continuación se muestra un ejemplo para n = 2. Es decir, en cada
sub-intervalo [yi , yi+1 ] se interpola con una función lineal (polinomio de grado n − 1 = 1)
que coincide en los extremos de ese sub-intervalo.
59
3.5. Ejercicios
Para tener en cuenta: en OCTAVE, un vector p de n componentes representa el
siguiente polinomio de grado n − 1:
>> p = [6 -3 0 1];
>> integrarpol(p)
ans =
1.5000 -1.0000 0 1.0000 0
>> p = [6 -3 0 1];
>> integrarpol(p,-3,5)
ans =
672
60
(a) Calcular el polinomio p(x) de grado 3 cuya gráfica pasa por esos puntos.
(a) Calcular el polinomio p(x) de grado 3 cuya gráfica pasa por esos puntos.
3.5. Considerar los siguientes números de la oficina de Censos de Estados Unidos, acerca
de la población de Estados Unidos:
Año Población
1900 75.994.575
1910 91.972.266
1920 105.710.620
1930 122.775.046
1940 131.669.275
1950 150.697.361
1960 179.323.175
1970 203.235.298
Graficar el polinomio usando un paso de 1 año entre 1900 y 1970. Repetir los gráficos
entre 1900 y 1980. El número relevado para 1980 es 226.547.082. ¿Qué opina sobre
la predicción dada por el polinomio interpolador?
61
x 0 1 2 3 4 5 6 7 8 9
y −0.1 1.1 1.9 3.2 3.8 5 6 7.3 8.1 8.9
x -1 0 1 3 6
y 6.1 2.8 2.2 6 26.9
x −1 0 2 3
y 0.3 −0.2 7.3 23.3
con funciones del tipo indicado en cada inciso. Graficar sobre el intervalo [−2, 4] las
curvas y los puntos.
1. y = a2x + b3x
2. y = a2x + b3x + c
3.10. Aproximar los datos de la tabla siguiente con un modelo de la forma f (x) ∼
=
ax2 +bx+c
−e , en el sentido de cuadrados mı́nimos para la función ln(f (x)):
x −1 0 1 2
y −1.1 −0.4 −0.9 −2.7
62
(a) Determinar la función lineal que mejor aproxima en el sentido de cuadrados mı́nimos
los datos dados (modelo lineal).
(c) Predecir cuál será la cantidad de mosquitos al cabo de 10 semanas según los diferentes
modelos, el lineal y el cuadrático.
(d) Si la medición a las 10 semanas es de 14900 mosquitos, ¿qué modelo le parece que
es más apropiado? ¿El lineal o el cuadrático?
63
Capı́tulo 4
El objetivo de este capı́tulo es hallar fórmulas que permitan calcular integrales (apro-
ximadas) pero de cualquier función. Cuando decimos aproximadas queremos decir con
tanta precisión como se desee.
Queremos calcular por ejemplo
Z A Z 5
−x2 t3
e dx o t
dt
0 6 e −1
Rara vez el integrando de esta expresión se podrá integrar analı́ticamente, y son muy
útiles los métodos numéricos. Estos no dan una fórmula cerrada para el resultado, pero
permiten calcular la integral con tanta precisión como se desee.
64
Esta fórmula se conoce como Regla del Rectángulo, porque estamos calculando el área del
rectángulo con base (b − a) y altura f ( a+b
2
).
Para n = 2 podemos escribir el polinomio p2 que interpola a f en x = a y en x = b
de la siguiente manera
b−x x−a
p2 (x) = f (a) + f (b).
b−a b−a
Luego
Z b Z b
1 1
f (x) dx ≈ p2 (x) dx = · · · = (b − a) f (a) + f (b) .
a a 2 2
| {z }
Q2 (f ;a,b)
Esta fórmula se conoce como Regla del Trapecio, porque estamos calculando el área del
trapecio con bases f (a), f (b) y altura (b − a).
Para n = 3 podemos escribir el polinomio p3 que interpola a f en x = a, x = m = a+b
2
y en x = b de la siguiente manera
Luego
Z b Z b
1 4 1
f (x) dx ≈ p3 (x) dx = · · · = (b − a) f (a) + f (m) + f (b) .
a a 6 6 6
| {z }
Q3 (f ;a,b)
(1) Las fórmulas de integración numérica consisten en evaluar f en distintos puntos del
intervalo, multiplicar esos valores por otros números, sumar, y finalmente multiplicar
por la longitud del intervalo.
65
Para hallar los wi correspondientes a la fórmula de Newton-Côtes de n puntos, consi-
deramos el intervalo [0, 1] y recordamos que deberı́a cumplirse que
Z b n
X
f (x) dx = Qn (f, a, b) = wi f (xi )
a i=1
w1 + w2 + · · · + wn = 1
1
w1 x1 + w2 x2 + · · · + wn xn =
2
1
w1 x21 + w2 x22 + · · · + wn x2n =
3
..
.
1
w1 xn−1
1 + w2 x2n−1 + · · · + wn xnn−1 = .
n
En forma matricial
1 1 ...
1
x1 w1
x2 ...
1
xn
w2 1
x21 x22
2... x2n
.. = ..
..
. .
. 1
n−1 n−1 n−1 wn n
x1 x2 . . . xn
Por lo tanto, los pesos de la fórmula de Newton-Côtes se encuentran resolviendo ese
sistema. Hay otras maneras de hallar los pesos, que son más precisas si n es grande, pero
para n ≤ 10 (que es lo que utilizaremos nosotros), esta manera es práctica.
El siguiente código calcula las abscisas xi y los pesos wi para el intervalo [0, 1] y un
entero positivo n.
66
Función pesosNC.m.
function [w,x] = pesosNC(n)
% function [w,x] = pesosNC(n)
% se calculan las abscisas y los pesos
% de la formula de Newton-Cotes de n puntos
x = linspace(0,1,n);
A = ones(n,n);
for i=2:n
A(i,:) = A(i-1,:) .* x;
end
b = 1./(1:n);
b = b’;
w = A\b;
El siguiente código calcula la integral aproximada de una función f utilizando la
fórmula de Newton-Côtes de n puntos. Las abscisas se calculan tomando n puntos equies-
paciados en el intervalo [a, b] (usando linspace).
Función integralNC.m.
function Q = integralNC(f,a,b,n)
% function Q = integralNC(f,a,b,n)
% se calcula la integral de f sobre [a,b]
% aproximadamente, usando la formula
% de Newton-Cotes de n puntos
[w,t] = pesosNC(n);
x = a+t*(b-a);
fx = f(x);
fx = fx(:);
w = w(:);
Q = (b-a)*sum(w .* fx);
67
con pn el polinomio de grado ≤ n − 1 que interpola a f en n puntos. Este polinomio
satisface, según el Corolario 41
Mn |b − a|n
máx |f (x) − pn (x)| ≤ ,
x∈[a,b] n (n − 1)n
con Mn = máxξ∈[a,b] |f n (ξ)|.
Por lo tanto
Z b Z b Z b
f (x) dx − Q n (f ; a, b) =
f (x) dx − p n (x) dx
a
Za b a
Z b
= f (x) − pn (x) dx ≤ |f (x) − pn (x)| dx
a a
n
Mn |b − a| Mn |b − a|n+1
≤ |b − a| =
n (n − 1)n n (n − 1)n
Hemos demostrado entonces el siguiente teorema:
Teorema 44. Si f es una función C n en [a, b], entonces
Z b
Mn |b − a|n+1
f (x) dx − Q n (f ; a, b)≤ ,
a
n(n − 1)n
con Mn = máxξ∈[a,b] |f n (ξ)|.
Esta cota del error entre la integral exacta y la aproximada tiene el mismo problema
que la cota del error de interpolación: Hay muchas funciones para las que ese término no
tiende a cero (cuando Mn → ∞ rápidamente.) Este problema se resuelve utilizando las
fórmulas de Newton-Côtes compuestas.
Luego, aproximamos cada integral con una fórmula de Newton-Côtes con n fijo:
Z yi+1
f (x) dx ≈ Qn (f, yi , yi+1 )
yi
y resulta
Z b L
X
f (x) dx ≈ Qn (f, yi , yi+1 ) .
a i=1
| {z }
QL
n (f ;a,b)
Para calcular una integral con este método, podemos usar el siguiente código.
68
Función NCcompuesta.m.
function Q = NCcompuesta(f,a,b,L,n)
% function Q = NCcompuesta(f,a,b,L,n)
% aproxima la integral de f sobre [a,b]
% utilizando la formula de Newton-Cotes compuesta
% de n puntos, subdividiendo en L subintervalos
y = linspace(a,b,L+1);
Q = 0;
for i=1:L
Q = Q + integralNC(f,y(i),y(i+1),n);
end
Esta implementación es muy ineficiente, pues cada vez que se llama a integralNC se
calculan los pesos de Newton-Côtes, y estos son los mismos en todos los sub-intervalos.
La siguiente, es una implementación un poco más eficiente.
Función intNCcompuesta.m.
function Q = intNCcompuesta(f,a,b,L,n)
% function Q = intNCcompuesta(f,a,b,L,n)
% aproxima la integral de f sobre [a,b]
% utilizando la formula de Newton-Cotes compuesta
% de n puntos, subdividiendo en L subintervalos
y = linspace(a,b,L+1);
h = y(2)-y(1);
Q = 0;
for i=1:L
x = y(i)+t*h;
fx = f(x); fx = fx(:);
Q = Q + h*sum(w .* fx);
end
69
4.4. Análisis del error para las fórmulas de Newton-
Côtes compuestas
Observemos lo siguiente:
Z b L Z yi+1 L
L
X X
f (x) dx − Q (f, a, b)= f (x) dx − Qn (f, yi , yi+1 )
n
a i=1 yi i=1
L Z yi+1
X
≤
f (x) dx − Q n (f, y i , yi+1 )
i=1 yi
n+1 n+1
L |b−a| |b−a|
X Mn L Mn Mn |b − a|n+1 1
L
≤ = L = .
i=1
n (n − 1)n n (n − 1)n n (n − 1)n Ln
Aquı́ hemos usado, en cada sub-intervalo [yi , yi+1 ] la estimación del error dada por el
Teorema 44. Hemos entonces demostrado el siguiente teorema.
Teorema 46. Si f es una función C n en [a, b], entonces
Z b
Mn |b − a|n+1 1
L
f (x) dx − Qn (f ; a, b) ≤ ,
a n(n − 1)n Ln
con Mn = máxξ∈[a,b] |f n (ξ)|.
Observación 47. Observemos que si f ∈ C n , entonces Mn es un número finito, y el lado
derecho de la fórmula anterior tiende a cero cuando L → ∞. Es decir, el error tiende a
cero a medida que tomamos L cada vez más grande. Además, cuando mayor sea n, más
rápido tenderá a cero el error.
Pero ¡cuidado! El error no tiende a cero si dejamos L fijo y hacemos tender n a infinito.
¿Por qué? Pensar.
4.5. Ejercicios
4.1. Utilizar la función pesosnc para determinar los pesos de la fórmula de Newton-Côtes
y los puntos xi para n = 2, n = 3, . . . , n = 6.
4.2. Consideremos Qn (f, 0, 5) la fórmula de Newton-Cotes de n puntos en el intervalo
R 5 5]. Para las siguientes funciones f : R → R llenar el cuadro y verificar si el error entre
[0,
0
f (x) dx y Qn (f ) tiende a cero a medida que n crece. Si no ocurre, pensar cuál puede
ser el motivo. También pensar por qué el error da tan grande en algunos casos y por qué
da cero en otros.
Z 5
n f (x) dx − Qn (f, 0, 5)
0
2
3
4
..
.
8
70
Z
1 1
(a) f (x) = , dx = tan−1 (x) + C
1 + x2 1 + x2
(f ) f (x) = |x|3/2 .
4.3. Consideremos QL2 (f, 0, 5) la fórmula del trapecio compuesta con parámetro de discre-
5
tización h = en el intervalo [0, 5]. Para cada una de las funciones del ejercicio anterior,
L
llenar un cuadro como el siguiente y verificar si el error tiende a cero y si se verifica que
el error para cada h es aproximadamente un cuarto del error para el h anterior. Si no
ocurre, pensar cuál puede ser el motivo.
5
Z
h
f (x) dx − QL2 (f, 0, 5)
0
1/2
1/4
1/8
1/16
..
.
1/1024
4.4. Repetir el ejercicio anterior con la fórmula de Simpson compuesta QL3 . ¿Cómo se
reduce el error a medida que h se divide por la mitad, para cada ejemplo?
Z 1
2
4.5. Se quiere calcular e−x dx utilizando la regla del trapecio compuesta, partiendo
−1
el intervalo [−1, 1] en L subintervalos. Hallar L de modo que el error sea menor que 10−3 .
Calcular la aproximación de la integral.
4.6. Programar la función I = trapcomp(x,y) que reciba como argumentos una lista de
valores de x y una lista de valores de y = f (x) y calcule la aproximación de la integral
de f en el intervalo que va desde x(1) hasta x(end) utilizando la regla del trapecio
compuesta (se supone que los xi están equiespaciados).
Utilizar la función para aproximar la integral de una función con estos valores:
71
x y
0 0
0.1571 0.1564
0.3142 0.3090
0.4712 0.4540
0.6283 0.5878
0.7854 0.7071
0.9425 0.8090
1.0996 0.8910
1.2566 0.9511
1.4137 0.9877
1.5708 1.0000
1.7279 0.9877
1.8850 0.9511
2.0420 0.8910
2.1991 0.8090
2.3562 0.7071
2.5133 0.5878
2.6704 0.4540
2.8274 0.3090
2.9845 0.1564
3.1416 0.0000
4.7. Programar la función I = simpsoncomp(x,y) que reciba como argumentos una lista
de valores de x y una lista de valores de y = f (x) y calcule la aproximación de la integral
de f en el intervalo que va desde x(1) hasta x(end) utilizando la regla de Simpson
compuesta (se supone que los xi están equiespaciados y las listas tienen un número impar
de valores).
Utilizar la función para aproximar la integral de una función con los valores de la
tabla del ejercicio anterior.
4.8. Se tiene la siguiente tabla de datos de una función f :
x y
0.0 1.0000000
0.5 0.4444444
1.0 0.1666667
1.5 0.0816327
2.0 0.0476190
2.5 0.0310078
3.0 0.0217391
3.5 0.0160643
4.0 0.0123457
4.5 0.0097800
5.0 0.0079365
Se sabe que |f 0 (x)| ≤ 1.46, que |f 00 (x)| ≤ 10 y que |f 000 (x)| ≤ 52.
R5
Se quiere aproximar la integral 0 f (x) dx usando esa tabla de valores.
72
(a) Dar una cota para el error cometido por la fórmula del trapecio compuesta.
(b) Dar una cota para el error cometido por la fórmula del Simpson compuesta.
4.9. Se tiene la siguiente tabla de valores de f (x). Se sabe que |f 00 (x)| ≤ 12, |f 000 (x)| ≤ 12
y |f (IV ) |(x)| ≤ 4 para todo x ∈ [0, 1].
x y
0 1.0000
0.0500 0.9998
0.1000 0.9980
0.1500 0.9932
0.2000 0.9840
0.2500 0.9688
0.3000 0.9460
0.3500 0.9143
0.4000 0.8720
0.4500 0.8177
0.5000 0.7500
0.5500 0.6672
0.6000 0.5680
0.6500 0.4508
0.7000 0.3140
0.7500 0.1562
0.8000 -0.0240
0.8500 -0.2282
0.9000 -0.4580
0.9500 -0.7147
1.0000 -1.0000
R1
¿Qué fórmula da menor error entre 0 f (x) dx y Q(f, 0, 1)?
¿La del trapecio compuesta o la de Simpson compuesta?
4.10. Calcular, con un error absoluto menor a 10−3 , el área de las superficies de revolución
que se obtienen al girar la gráfica de y = f (x) alrededor del eje x, para x en el intervalo
indicado.
A modo de ejemplo, mostramos las lı́neas que deben escribirse en MATLAB para graficar
la superficie de revolución del ı́tem (a).
73
x = linspace(0,2,20);
phi = linspace(0,2*pi,30);
[x,phi] = meshgrid(x,phi);
f = @(x) 2+cos(pi*x);
y = cos(phi).*f(x);
z = sin(phi).*f(x);
surf(x,y,z)
74
Capı́tulo 5
5.1. Introducción
El objetivo en un problema a valores iniciales (PVI) es encontrar una función y(t)
dado su valor en un instante inicial t = 0 y una receta f (t, y) para su pendiente:
(
y 0 = f (t, y), T0 ≤ t ≤ Tf ,
(PVI)
y(T0 ) = y0 .
En algunas aplicaciones se desea graficar una aproximación a y(t) sobre un intervalo de
interés [T0 , Tf ] para descubrir propiedades cualitativas de la solución. En otras, puede ser
importante conocer una estimación precisa del valor y(t) para un tiempo t = Tf prefijado.
Los datos conocidos de (PVI) son: el valor inicial y0 y la función pendiente f (t, y) que
nos dice cuál tiene que ser la pendiente de la solución y(t) en cada instante de tiempo,
dependiendo de cuánto vale y. Más precisamente:
Se llama solución de (PVI) a una función y(t) definida sobre el intervalo [T0 , Tf ] que
satisface:
y(T0 ) = y0 , y 0 (t) = f (t, y(t)), para todo t ∈ [T0 , Tf ].
Consideremos por ejemplo el PVI
(
y 0 = −2y, 0 ≤ t < ∞,
y(0) = 3.
En este caso, el valor inicial es 3 y la función pendiente es f (t, y) = −2y (es una función
de dos variables que en realidad depende de una sola). La ecuación diferencial y 0 = −2y
tiene una familia de soluciones, son todas las funciones de la forma y(t) = ce−2t , con c
cualquier constante. La condición inicial y(0) = 3 nos dice que queremos hallar la única
función de esa familia que en t = 0 vale 3. Utilizamos la condición inicial para despejar
c y nos da que c = 3. Luego, la solución del PVI es y(t) = 3e−2t .
75
Si f es una función continua de sus dos variables (y y t) y además ∂f
∂y
existe y está
acotada, entonces el problema a valores iniciales (PVI) tiene una única solución. Más
aún, si L ≤ ∂f
∂y
≤ U para todo (y, t), resulta lo siguiente:
|y0 − ỹ0 |eL(t−T0 ) ≤ |y(t) − ỹ(t)| ≤ |y0 − ỹ0 |eU (t−T0 ) , para t ∈ [T0 , Tf ].
En palabras, el resultado anterior nos dice que las soluciones del mismo problema, con
diferentes condiciones iniciales, se alejan al menos tan rápido como eLt y no más rápido
que eU t por la diferencia entre las condiciones iniciales.
En el ejemplo anterior, f (t, y) = −2y, luego ∂f
∂y
= −2, y por lo tanto L = U = −2. En
este caso, tenemos que soluciones del mismo PVI con dos condiciones iniciales diferentes
se acercan como e−2t por el error inicial. Más precisamente, observemos las soluciones
con y0 = 3 y con ỹ0 = 3.1. En este caso,
por lo que
|y(t) − ỹ(t)| = |3e−2t − 3.1e−2t | = |3 − 3.1|e−2t = 0.1e−2t .
76
En palabras, conocido (o aproximado) el valor de y(t), el valor de y(t + h) es aproxima-
damente y(t) + hf (t, y(t)).
Definir ti = (i − 1) h, i = 1, 2, . . . L + 1.
Es decir:
t1 = T0 Y1 = y0
t2 = T0 + h Y2 = Y1 + h f (t1 , Y1 )
t3 = T0 + 2h Y3 = Y2 + h f (t2 , Y2 )
t4 = T0 + 3h Y4 = Y3 + h f (t3 , Y3 )
.. ..
. .
tL+1 = T0 + Lh = Tf YL+1 = YL + h f (tL , YL )
Figura 5.1: Solución de (PVI) para f (t, y) = −2y, y0 = 3 sobre el intervalo [0, 4] y compa-
ración con el método de Euler. La lı́nea curva representa la solución exacta y(t) = 3e−2t .
Los asteriscos son los puntos (ti , Yi ), y los cı́rculos los puntos (ti , y(ti )). A la izquierda
la solución aproximada se obtuvo con L = 20 (h = 0.2) y el error global es 0.268. A la
derecha, la solución aproximada se obtuvo con L = 40 (h = 0.1) y el error global es 0.121
77
Función euler1.m.
function [t,y] = euler1(f,inter,y0,L)
% function [t,y] = euler1(f,[T0 TF],y0,L)
% Metodo de Euler para resolver
% y’ = f(t,y) en [t0,TF]
% y(t0) = y0
% Usando L pasos
t = linspace(inter(1),inter(2),L+1)’;
h = (inter(2)-inter(1))/L;
y(1) = y0;
for j = 1:L
y(j+1) = y(j) + h*f(t(j),y(j));
end
Es claro que en cada paso se comete un error, pequeño si h es pequeño, pero estos
errores de alguna manera se acumulan, y para llegar a Tf la cantidad de pasos será mayor
cuanto menor sea h.
Para entender un poco mejor el error cometido, necesitamos definir algunas cantidades:
Error global: Se llama error global al máximo error entre y(tn ) y Yn , es decir,
|yj (t) − yj+1 (t)| ≤ |yj (tj+1 ) − yj+1 (tj+1 )| = |yj (tj+1 ) − Yj+1 | = εj .
78
Entonces
gn = |Yn − y(tn )| = |Yn − y1 (tn )|
≤ |Yn − yn−1 (tn )| + |yn−1 (tn ) − yn−2 (tn )| + · · · + |y2 (tn ) − y1 (tn )|
n−1
X
≤ εn−1 + εn−2 + · · · + ε1 = εj .
j=1
Hasta aquı́, el análisis que hemos hecho es independiente del método, si es Euler o
algún otro. Veamos cómo se puede acotar εj para el método de Euler.
Una herramienta muy utilizada para estimar y acotar el error en los métodos numéricos
es el Teorema de Taylor, que enunciamos a continuación:
Teorema 48 (Teorema de Taylor). Sea y ∈ C n+1 (I), donde I es un intervalo que contiene
a t0 . Entonces, para cada t ∈ I existe x entre t0 y t tal que
(t − t0 )2 00 (t − t0 )n (n) (t − t0 )n+1 (n+1)
y(t) = y(t0 ) + (t − t0 )y 0 (t0 ) + y (t0 ) + · · · + y (t0 ) + y (x).
2! n! (n + 1)!
Tf −T0
Usando el teorema de Taylor y considerando que h = L
y que tj+1 = tj + h,
tenemos que
h2 00 h2
yj (tj+1 ) = yj (tj ) + hyj0 (tj ) + yj (ξj ) = Yj + hf (tj , Yj ) + yj00 (xj )
2 | {z } 2
Yj+1
Concluimos entonces que si todas las soluciones de la EDO tienen derivada segunda
acotada uniformemente por M , el método de Euler produce una solución con máximo
M (Tf − T0 )
error acotado por h = Const. × h. Se dice que el método de Euler es de orden
2
uno (la potencia de h). La demostración que hemos visto es muy simple, y analizando el
∂f
caso ≤ 0, permite entender la idea general que si un método produce en un paso un
∂y
error de orden hk+1 entonces el error global será de orden hk .
79
5.4. PVI para sistemas de ecuaciones y ecuaciones
de orden superior
Consideremos ahora un sistema de ecuaciones diferenciales ordinarias. Por ejemplo,
un modelo de predador presa:
(
x01 = 0.5x1 + 0.2x1 x2 x1 (0) = 500
0
x2 = 0.8x2 − 0.8x1 x2 x2 (0) = 10000
g
ϕ00 (t) = − sen ϕ(t)
`
donde g denota la constante de aceleración gravita-
cional. Esta es una ecuación diferencial ordinaria de
orden dos, pues aparecen derivadas segundas de la
incógnita ϕ.
80
Este problema resulta un PVI cuando se complementa con dos condiciones iniciales:
ϕ(0) = ϕ0 , ϕ0 (0) = v0 ,
y1 = ϕ,
y2 = ϕ 0 ,
de manera que
y10 = ϕ0 = y2
g g
y20 = ϕ00 = − sen ϕ = − sen y1 ,
` `
y
y1 (0) = ϕ(0) = ϕ0 , y2 (0) = ϕ0 (0) = v0 .
Ahora reescribimos las ecuaciones en términos de y1 , y2 sin usar ϕ, ϕ0 .
y10 = y2
(
y1 (0) = ϕ0 ,
0 g
y2 = − sen y1 , y2 (0) = v0 .
L
y
Si llamamos y = 1 resulta
y2
0 y2 ϕ0
y = y(0) = .
− Lg sen y1 v0
[T0 TF]: Es un vector de dos componentes que indican el intervalo donde se desea
conocer la solución.
81
Salidas:
t: Debe ser un vector columna de L+1 componentes.
y: Debe ser una matriz de L+1 filas y n columnas. La i-ésima fila debe contener el
estado y a tiempo t(i).
Función euler.m.
function [t,y] = euler(f,inter,y0,L)
% function [t,y] = euler(f,[t0 TF],y0,L)
% Metodo de Euler para resolver
% y’ = f(t,y) en [t0,TF]
% y(t0) = y0
% Usando L pasos
% y0 puede ser vectorial o escalar
t = linspace(inter(1),inter(2),L+1)’;
h = (inter(2)-inter(1))/L;
y(1,:) = y0’;
for i = 2:L+1
y(i,:) = y(i-1,:) + h*f(t(i-1),y(i-1,:)’)’;
end
82
Por el teorema de Taylor,
h2 00 h3
y(tn+1 ) = y(tn ) + hy 0 (tn ) + y (tn ) + y 000 (x),
2 3!
pero
h2 n h3
y(tn+1 ) = Yn + hf n + ft + fyn f n + y 000 (x).
2 3!
Por otro lado
Yn+1 = Yn + ak1 + bk2
con
k1 = hf n
k2 = hf (tn + αh, Yn + αk1 )
= h f (tn ) + αhft (tn , Yn ) + αk1 fy (tn , Yn ) + O(h2 )
Por lo tanto
Yn+1 = Yn + (a + b)hf n + bh2 αftn + αf n fyn + O(h3 ).
De esta manera,
1
Yn+1 − y(tn+1 ) = Yn − Yn + [(a + b) − 1]hf n + bα(ftn + f n fyn ) − (ftn + f fyn ) h2 + O(h3 ),
2
y será εn = |Yn+1 − y(tn+1 )| = O(h3 ) si
a+b=1
1
bα =
2
Elegimos la solución simétrica
1 1
a= , b= , α = 1,
2 2
y obtenemos el método de Runge-Kuta de orden 2:
83
Dados tn e Yn , calcular tn+1 e Yn+1 a partir de las siguientes fórmulas:
k1 = hf (tn , Yn )
k2 = hf (tn + h, Yn + k1 )
(RK2) 1 1
Yn+1 = Yn + k1 + k2
2 2
tn+1 = tn + h.
Teorema: Si todas las soluciones tienen derivada de orden tres acotada uni-
formemente por una constante M3 , entonces el método de Runge-Kutta de
orden 2 produce una sucesión (tn , Yn ) que cumple
T −T
donde h = f L 0 y C depende de M3 y de la longitud del intervalo [T0 , Tf ], es
decir (Tf − T0 ).
De manera análoga, pero con cuentas más engorrosas puede deducirse el método de
Runge-Kutta de orden cuatro:
Cabe aclarar aquı́ que el método de Runge-Kutta de orden cuatro convergerá con ese
orden cuando las derivadas de hasta orden cinco de y (y cuatro de f ) estén acotadas.
Más precisamente Para este método vale el siguiente resultado:
T −T
donde h = f L 0 y C depende de M4 y de la longitud del intervalo [T0 , Tf ], es
decir (Tf − T0 ).
84
La pregunta que cabe hacerse es la siguiente: ¿Cuál es la ganancia de usar un método
de orden superior si es computacionalmente más costoso en cada paso de tiempo? Para
ver la ganancia, consideremos el caso en que las cotas del error de Euler, RK2 y RK4 son
h, h2 y h4 respectivamente, es decir, C = 1 en todos los casos. Recordemos que la parte
más costosa desde el punto de vista computacional es la evaluación de f (t, y), por lo que
interpretamos que Euler, RK2 y RK4 cuestan 1, 2 y 4 “operaciones” por paso de tiempo.
En la siguiente tabla mostramos el costo y el error de cada método para diferentes L
pensando en el intervalo [0, 1].
k1 = hf (tn , Yn ),
j−1
X
kj = hf (tn + αj h, Yn + βji ki ), j = 2, 3, . . . , 6,
i=1
1
Los valores pueden encontrarse en las tablas 6.25 y 6.26 de la página 430 del libro de Atkinson
85
con h = hn+1 = tn+1 − tn , y
5
X
Yn+1 = Yn + γj kj ,
j=1
6
X
Ŷn+1 = Yn + γ̂j kj ,
j=1
donde (
yn0 = f (t, yn ),
yn (tn ) = Yn .
Luego
La idea del método adaptativo es la siguiente: Fijar los parámetros 0 < A << 1,
0 < P < 1 < Q, y una tolerancia tol; en cada paso hacer lo siguiente:
86
Observación 51. Cuando el algoritmo termina, digamos, en N pasos, resulta que
N N N
X X X hn tol
εn ≈ Tn ≤ = tol,
n=0 n=0 n=0
b−a
y por lo tanto,
máx |y(tn ) − Yn | ≤ Ctol,
0≤n≤N +1
∂f
donde C depende del máximo valor de y de la longitud del intervalo.
∂y
El paso marcado con (*) sirve para agrandar el paso de discretización hn cuando el
error es mucho más pequeño que lo deseado. Esto permite ahorrar operaciones.
El paso marcado con (**) permite reducir el paso de discretización hn cuando el
error es más grande que lo permitido.
Observación 52. La idea de este método RKF45 es la base de la función ode45 de
MATLAB, que tiene además otros algoritmos adaptativos para EDOs: ode23, ode15s,
ode23s. Los que terminan en s están diseñados para problemas stiff, que usan métodos
∂f
implı́citos y son especialmente aptos para problemas donde tiene autovalores muy
∂y
negativos.
cuya solución exacta es y(t) = e−100t . Este es un problema estable dado que la solución
tiende a cero cuando t → ∞ y más aún, si tomamos soluciones con y(0) = y0 y ỹ(0) = ỹ0 ,
entonces se cumple
|y(t) − ỹ(t)| ≤ e−100t |y0 − ỹ0 |.
Es decir, la diferencia existente a tiempo t = 0 disminuye (¡exponencialmente!) cuando t
crece.
Consideremos ahora su resolución con el método de Euler, con h = 0.1. Resulta:
Y1 =1
Y2 = Y1 + hf (t1 , Y1 ) = Y1 + h(−100Y1 ) = 1 + 0.1(−100 × 1) = 1 − 10 = −9
Y3 = Y2 + hf (t2 , Y2 ) = Y2 + h(−100Y2 ) = −9 + 0.1(−100 × (−9)) = −9 + 90 = 81
Y4 = Y3 + hf (t3 , Y3 ) = Y3 + h(−100Y3 ) = 81 + 0.1(−100 × 81) = 81 − 810 = −729
Y5 = · · · = 6561
Y6 = · · · = −59049
..
.
¡Nada que ver con la solución!
87
Probemos con h más pequeño, por ejemplo h = 0.05. Resulta
Y1 =1
Y2 = Y1 + hf (t1 , Y1 ) = Y1 + h(−100Y1 ) = 1 + 0.05(−100 × 1) = 1 − 5 = −4
Y3 = Y2 + hf (t2 , Y2 ) = Y2 + h(−100Y2 ) = −4 + 0.05(−100 × (−4)) = −4 + 20 = 16
Y4 = Y3 + hf (t3 , Y3 ) = Y3 + h(−100Y3 ) = 16 + 0.05(−100 × 16) = 16 − 80 = −64
Y5 = · · · = 256
Y6 = · · · = −1024
..
.
¡Nada que ver con la solución!
Y1 =1
Y2 = Y1 + hf (t1 , Y1 ) = Y1 + h(−100Y1 ) = 1 + 0.02(−100 × 1) = 1 − 2 = −1
Y3 = Y2 + hf (t2 , Y2 ) = Y2 + h(−100Y2 ) = −1 + 0.02(−100 × (−1)) = −1 + 2 = 1
Y4 = Y3 + hf (t3 , Y3 ) = Y3 + h(−100Y3 ) = 1 + 0.02(−100 × 1) = 1 − 2 = −1
Y5 = ··· = 1
Y6 = · · · = −1
..
.
¡Nada que ver con la solución!
Y1 =1
Y2 = Y1 + hf (t1 , Y1 ) = Y1 + h(−100Y1 ) = 1 + 0.01(−100 × 1) = 1 − 1 = 0
Y3 = Y2 + hf (t2 , Y2 ) = Y2 + h(−100Y2 ) = 0 + 0.02(−100 × 0) = −1 + 2 = 0
Y4 = Y3 + hf (t3 , Y3 ) = Y3 + h(−100Y3 ) = 0 + 0.02(−100 × 0) = 1 − 2 = 0
Y5 = ··· = 0
Y6 = ··· = 0
..
.
Esta ya se parece un poco...
Para h más pequeños ya se parece un poco más. Pero para h > 0.01 las soluciones
son realmente diferentes a la solución exacta y = e−100t . Es más, si h > 0.02, Yn → ∞
cuando n → ∞.
El concepto de estabilidad más estudiado en métodos numéricos para ecuaciones di-
ferenciales tiene que ver con el comportamiento del método al resolver el siguiente PVI
para λ ∈ C: (
y 0 = λy, t > 0,
y(0) = 1,
88
cuya solución exacta para λ = λR + iλI es
λt λR t
y(t) = e = e cos λI t + i sin λI t .
Y1 = 1,
Yn+1 = Yn + hλYn = (1 + hλ)Yn = (1 + hλ)n+1 ,
por lo que la sucesión {Yn }∞n=1 será acotada si y solo si |1 + hλ| ≤ 1. Por lo tanto, la
región de estabilidad para el método de Euler es
Observación 54.
Si λ < 0, h debe satisfacer h ≤ −2/λ, ası́ que debe ser muy pequeño cuando |λ| es
grande.
89
Cuando se quieren diseñar métodos más estables, se entiende que se desean métodos
con una mayor región de estabilidad.
∂f
Cuando un problema y 0 = f (t, y) tiene con autovalores de parte real muy
∂y
negativa (que exige h muy pequeño) se dice que la EDO es stiff (rı́gida). Esta
situación es muy común en problemas que involucran reacciones quı́micas.
k1 = hλYn ,
k2 = hλ(Yn + k1 ) = hλ(Yn + hλYn ) = [hλ + (hλ)2 ]Yn ,
hλ hλ (hλ)2
k1 + k2
Yn+1 = Yn + = Yn 1 + + +
2 2 2 2
2
(hλ)
= Yn 1 + hλ + .
2
2 n
Luego, Yn+1 = 1 + hλ + (hλ)
2
y el método es estable si y solo si
z2
hλ ∈ RRK2 := {z ∈ C : |1 + z + | ≤ 1}.
2
La región R de estabilidad de RK2 puede verse en la Figura 5.2.
90
Region de estabilidad Euler Region de estabilidad RK4
Region de estabilidad RK2
3 3 3
2 2 2
1 1 1
0 0 0
−1 −1 −1
−2 −2 −2
−3 −3 −3
−4 −3 −2 −1 0 1 −4 −3 −2 −1 0 1 −4 −3 −2 −1 0 1
Figura 5.2: Regiones de estabilidad para métodos de Runge-Kutta. Euler o RK1 (izquier-
da), RK2 (medio), RK4 (derecha).
Y = Yn + hf (tn+1 , Y ).
Este método implica resolver, en cada paso, una ecuación o sistema de ecuaciones que
en general será no lineal. En cada paso se puede resolver con un método tipo Newton
(por ejemplo). Este método tiene un costo adicional de resolver un sistema no lineal en
cada paso de tiempo, pero permite tomar ∂f h mucho más grande que el método de Euler
∂f
explı́cito, en especial cuando < 0 y es grande (ecuaciones stiff). Veamos la región
∂y ∂y
de estabilidad de Retro Euler (RE). Si aplicamos la definición Yn+1 = Yn + hf (tn+1 , Yn+1 )
al problema modelo y 0 = λy obtenemos que
1 1 n+1
Yn+1 = Yn + hλYn+1 ⇒ Yn+1 = Yn = Y0 .
1 − hλ 1 − hλ
Por lo tanto, el método resultará estable cuando hλ pertenezca a la región de estabilidad
n 1 o
RRE = z ∈ C : ≤ 1 = {z ∈ C : |1 − z| ≥ 1} = B(1, 1)c ,
1−z
91
que contiene a todo el semiplano complejo {z ∈ C : Re z ≤ 0}. En este caso se dice que
el método es incondicionalmente estable, pues h no debe cumplir ninguna condición para
que el método sea estable, cuando Re λ ≤ 0.
5.10. Ejercicios
5.1. Verifique el funcionamiento del método de Euler (euler1.m), utilizándolo para re-
solver el siguiente PVI: (
y 0 = −y, 0 ≤ t ≤ 1,
y(0) = 1.
Este PVI tiene solución exacta y(t) = e−t . Llene la siguiente tabla:
L h error
10 0.1
20 0.05
40 0.025
80 0.0125
160 0.00625
320 0.003125
con el método de Euler. Utilice diferentes valores de L (por ejemplo: 10, 20, 40, 80, 160,
. . . ) y grafique las soluciones que se van obteniendo. Determine visualmente un valor de
L para el que la solución no cambia más. Finalmente calcule el máximo y el mı́nimo de
la solución en el intervalo [0, 4].
5.3. Utilizando el valor de L del ejercicio anterior, calcule las soluciones de la misma
ecuación con condición inicial y(0) = 1, y(0) = 0 y y(0) = −1. Grafı́quelas todas sobre el
mismo par de ejes y diga qué ocurre a medida que el tiempo t avanza. ¿Se acercan o se
alejan las trayectorias? ¿Cómo se relaciona esto con ∂f∂y
?
5.4. Utilizando L = 10000, resuelva el PVI
(
y 0 = y(1 − y)(1 + y), 0 ≤ t ≤ 2,
y(0) = y0 ,
con y0 = −1.4, −1.3, −1.2, . . . , 1.3, 1.4. Superponga todas las soluciones en un gráfico y
explique el comportamiento cualitativo de las soluciones.
92
5.5. Consideremos el siguiente PVI correspondiente a un sistema de EDOs:
( (
x01 (t) = −x2 (t) x1 (0) = 1
0 0 ≤ t ≤ 2π, .
x2 (t) = x1 (t) x2 (0) = 0
La solución exacta es x1 (t) = cos(t), x2 (t) = sen(t). Realice una tabla similar a la del
Ejercicio 1 y verifique si el error tiende a cero con orden O(h), es decir, verifique si el
error se divide aproximadamente por la mitad cada vez que h se divide por la mitad.
5.6. La trayectoria de una partı́cula que se mueve en el plano está dada por la curva
(x1 (t), x2 (t)), donde las funciones x1 y x2 son la solución del siguiente sistema de ecua-
ciones diferenciales: (
x01 (t) = −tx2 (t)
x02 (t) = tx1 (t) − tx2 (t)
Resuelva este sistema en el intervalo [0, 20] con el método de Euler utilizando paso h =
0.05 y grafique la trayectoria de la partı́cula, sabiendo que en tiempo t = 0 se encontraba
en el punto (1, −1).
5.7. Programe el método de Runge-Kutta de orden 2, el formato debe ser:
[t,y] = rk2(f, [t0 TF], y0, L)
Donde las entradas f, [T0 Tf], y0, L y las salidas [t,y] son como se describe en el
apunte, previo a la función euler.m.
Verificar su funcionamiento repitiendo el ejercicio 5.5, pero con el método de Runge-
Kutta de orden 2, en lugar del de Euler. Si el método está bien implementado, cada vez
que h se divide por la mitad, el error deberı́a dividirse por 4(= 22 ).
5.8. Repetir el ejercicio anterior para el método de Runge-Kutta de orden 4. Es decir,
programar la función rk4 con el formato:
[t,y] = rk4(f, [t0 TF], y0, L)
Verificar su funcionamiento repitiendo el ejercicio 5.5, pero con el método de Runge-
Kutta de orden 4, en lugar del de Euler. Si el método está bien implementado, cada vez
que h se divide por la mitad, el error deberı́a dividirse por 16(= 24 ).
5.9. El siguiente PVI de orden 3 tiene como solución exacta y(t) = cos(t):
y + 4y 00 + 5y 0 + 2y = −4 sen t − 2 cos t
(3)
y(0) = 1
y 0 (0) = 0
y 00 (0) = −1
Resuélvala sobre el intervalo [0, 6] con Euler, RK2 y RK4 y llene una tabla como la
siguiente con el error global:
h Euler RK2 RK4
1/10
1/20
1/40
1/80
1/160
1/320
93
Determine a partir de la tabla, el valor de h y el número de evaluaciones de f que fueron
necesarios, en cada método, para tener un error global:
(a) ¿Cuál variable representa al predador y cuál a la presa? Para determinar esto, puede
resolver el sistema en los siguientes casos:
(b) Calcule numéricamente una solución donde la población inicial de la presa es 1000
y la de predadores es 500. Dibujar la solución, graficando ambas poblaciones con el
tiempo, y describir el fenómeno representado (utilice un método que usted considere
apropiado, y un intervalo de tiempo que sea representativo de lo que sucede).
g
ϕ00 (t) +
sen ϕ(t) = 0,
`
donde g = 9.81g/m2 es la constante de aceleración gravitacio-
nal.
Supongamos (por simplicidad) que la longitud del brazo es igual a 9.81m, con lo que
obtenemos la ecuación
ϕ00 (t) + sen ϕ(t) = 0.
Resolver esta ecuación numéricamente (utilizando ode45) en el intervalo [0, 10], en los
siguientes casos, y explicar la situación fı́sica descripta (condiciones iniciales y evolución):
94
5.12. Un sistema quı́mico muy común en la industria de procesos es el Reactor en tan-
que continuamente agitado (CSTR por su nombre en inglés: Continuously Stirred Tank
Reactor ). En este trabajo estudiaremos un reactor en tanque enfriado2 . Se supone que
el recipiente está perfectamente mezclado, y que ocurre una reacción exotérmica e irre-
versible de primer grado A → B. Un diagrama esquemático del dispositivo y la chaqueta
enfriante se muestra en la figura.
95
Tj con la temperatura T (t) del tanque. El segundo término indica la influencia sobre la
temperatura del reactor causada por la reacción quı́mica. En esta ecuación, H es un
parámetro de calor de reacción, cp un término de capacidad calórica, ρ una densidad, U
un coeficiente de transferencia de calor, y A el área para el intercambio de calor (área de
contacto entre el refrigerante y el tanque).
Poniendo todo junto, el CSTR tiene tres señales de entrada:
CAf Concentración de A en el caudal de entrada [kgmol/m3 ]
Tf Temperatura de la alimentación [K]
Tj Temperatura del lı́quido refrigerante [K]
Dos señales de salida, o variables de estado:
y1 (t) = CA (t) Concentración de A en el tanque reactor [kgmol/m3 ]
y2 (t) = T (t) Temperatura del reactor [K]
Después de agrupar algunos de los parámetros originales, llegamos a ocho parámetros
diferentes:
F Tasa de flujo volumétrico [volumen/tiempo] [m3 /h] 1
3
V Volumen del reactor [m ] 1
k0 Factor no-térmico pre-exponencial [1/h] 35 × 106
E Energı́a de activación [kcal/kgmol] 11850
R Constante de gas de Bolzman [kcal/(kgmol K)] 1.98589
H Calor de la reacción [kcal/kgmol] -5960
3
HD = cp ρ Capacidad calórica por densidad [kcal/(m K)] 480
3
HA = U A Tasa total de transferencia de calor por área [kcal/(m K h)] 145
Con estas unidades, el tiempo t se mide en horas.
1. Escribir un script MATLAB para simular el reactor. Permitir que las constantes y
datos del problema sean fácilmente modificadas.
2. Simular el comportamiento durante un dı́a del reactor en los siguientes casos, y ela-
borar un pequeño informe con lo que se observa en las señales de salida, o variables
de estado (como para ser leı́do por un ingeniero).
Si le interesa leer un poco más sobre los temas de este capı́tulo, le recomendamos los
siguientes libros:
K.E. Atkinson, An introduction to Numerical Analysis, 2nd edition, John Wiley &
Sons, 1989.
97
Capı́tulo 6
Una ecuación en derivadas parciales (EDP) es una ecuación matemática que contiene
derivadas parciales, por ejemplo:
∂u ∂u
+3 = 0. (6.1)
∂t ∂x
Podrı́amos comenzar nuestro estudio determinando qué funciones u(x, t) satisfacen (6.1).
Sin embargo, preferimos comenzar investigando un problema fı́sico.
El estudio de las EDP domina diversas áreas de la ingenierı́a y de la fı́sica. He aquı́
una lista (no exhaustiva) de áreas altamente dependientes del estudio de EDP: acústica,
aerodinámica, elasticidad, electrodinámica, dinámica de fluidos, geofı́sica, transferencia
de calor, meteorologı́a, etc.
Con respecto al estudio de EDP seguiremos cierta filosofı́a de matemática aplicada en
la cual el análisis de un problema tendrá tres etapas:
1. Formulación
2. Resolución
3. Interpretación
En este curso estudiaremos las ecuaciones del flujo de calor que describen la trans-
ferencia de energı́a térmica. La energı́a térmica es causada por la agitación de materia
molecular. El movimiento de la energı́a térmica es causado por dos procesos básicos:
conducción y convección. La conducción resulta de las colisiones de moléculas vecinas
en las cuales la energı́a cinética de la vibración de una molécula es transmitida a la de
su vecina más cercana. De este modo la energı́a térmica se propaga por conducción aún
si las moléculas no se mueven de una región a otra. Además, si una molécula que está
vibrando se mueve de una región a otra, se lleva la energı́a térmica con ella. Este tipo de
movimiento de energı́a térmica se denomina convección. Para comenzar nuestro estudio
con problemas relativamente simples, estudiaremos el flujo de calor sólo en casos donde
la conducción de energı́a calórica es mucho más significante que su convección. Es decir,
estaremos pensando en flujo de calor en el caso de sólidos. Aunque la transferencia de
calor en fluidos (lı́quidos y gases) está también gobernada principalmente por conducción
si la velocidad del fluido es suficientemente pequeña.
98
6.1. Derivación de la Conducción del Calor en una
Barra Unidimensional
Densidad de Energı́a Térmica. Comenzamos considerando una barra de sección
transversal constante de área A, orientada en la dirección del eje x (desde x = 0 hasta
x = L) como se ve en la Figura 6.1. Introducimos la cantidad de energı́a térmica por
unidad de volumen como una variable desconocida (incógnita) y la llamamos densidad
de energı́a térmica:
e(x, t) ≡ densidad de energı́a térmica.
Supondremos de ahora en adelante que todas las cantidades térmicas son constantes
en cualquier sección transversal (esto quiere decir que consideramos la barra como uni-
dimensional). La forma más sencilla de lograr esto es aislar perfectamente la superficie
lateral de la barra. De este modo no hay energı́a pasando a través de la superficie lateral
de la barra. La dependencia de x y t corresponde al caso en que la distribución de
temperatura no es uniforme; la energı́a térmica varı́a de una sección transversal a otra.
Figura 6.1: Energı́a calórica fluyendo hacia dentro y hacia fuera de un segmento de la
barra.
99
térmicamente aislada. El proceso fundamental de flujo de calor se describe entonces por
la ecuación
razón de cambio Z b Z b
d d
de la energı́a = A e(x, t) dx = A e(x, t) dx.
dt a dt a
en el tiempo
Rb
Técnicamente, aquı́ aparece una derivada ordinaria d/dt, pues la expresión a
e(x, t) dx
depende sólo de t y no depende de x. Por otro lado, para a y b fijos,
Z b Z b
d ∂e
e(x, t) dx = (x, t) dx.
dt a a ∂t
Flujo de Calor. En una barra unidimensional como la que estamos considerando noso-
tros, la energı́a térmica fluye hacia la derecha o hacia la izquierda. Introducimos ahora la
variable φ(x, t) que llamaremos flujo de calor y que contiene la cantidad de energı́a térmi-
ca que fluye hacia la derecha por unidad de tiempo por unidad de área. Si φ(x, t) < 0,
entonces la energı́a está fluyendo hacia la izquierda.
La energı́a calórica que fluye por unidad de tiempo hacia adentro de la sección entre
x = a y x = b está dada por
Ahora bien, como φ(x, t) es el flujo por unidad de área y por lo tanto debe multiplicarse
por el área para saber el flujo total a través de una sección, tenemos que
Luego
100
Fuentes de Calor. Nuestro modelo también nos permitirá considerar la existencia de
fuentes internas de energı́a térmica, y utilizaremos la variable Q(x, t) para denotar la
energı́a calórica generada por unidad de tiempo por unidad de volumen. Esta energı́a
generada puede deberse —por ejemplo— a reacciones quı́micas o a calentamiento por
medio de electricidad.
Como Q(x, t) representa la energı́a calórica generada por unidad de tiempo por unidad
de volumen, para obtener la energı́a generada por unidad de tiempo en la sección que va
desde x = a hasta x = b debemos integrar, y obtenemos que
Z b
energı́a generada dentro ∼
= A Q(x, t) dx.
de la secciónpor unidad de tiempo a
Cancelando A, resulta
Z b Z b
∂e
(x, t) dx = φ(a, t) − φ(b, t) + Q(x, t) dx. (6.3)
a ∂t a
Como todas las integrales están escritas sobre el mismo intervalo [a, b] resulta que
Z b
∂e ∂φ
(x, t) + (x, t) − Q(x, t) dx = 0.
a ∂t ∂x
Esta integral debe ser cero para toda sección a < x < b, y por lo tanto, el integrando
debe ser cero. Ası́ obtenemos la siguiente ecuación de balance:
∂e ∂φ
(x, t) = − (x, t) + Q(x, t). (6.4)
∂t ∂x
101
la temperatura —que denotaremos con u(x, t)— en términos de la densidad de energı́a
térmica e(x, t). Para poder relacionar la energı́a térmica con la temperatura necesitamos
introducir el concepto de calor especı́fico o capacidad calórica:
Ley de Fourier. Usualmente, la ecuación (6.6) se entiende como una ecuación con dos
incógnitas: la temperatura u(x, t) y el flujo de calor (flujo por unidad de área por unidad
de tiempo) φ(x, t). Surge ahora la pregunta: ¿Cómo y por qué fluye la energı́a calórica?
En otras palabras, ¿Cuál es la relación entre el flujo de energı́a calórica y el campo de
temperaturas?
102
Primero resumimos ciertas propiedades cualitativas del flujo de calor con las que ya
estamos familiarizados:
4. El flujo de energı́a calórica será diferente para diferentes materiales, aun para la misma
diferencia de temperatura.
En general pensamos que Q(x, t) es un dato dado, y resulta ser la temperatura u(x, t)
la única función incógnita. Los coeficientes térmicos c, ρ, K0 dependen del material y
103
pueden ser funciones de x. En el caso especial de una barra uniforme en que c, ρ, K0 son
todos constantes, la ecuación diferencial (6.8) se transforma en
∂u ∂ 2u
cρ (x, t) = K0 2 (x, t) + Q(x, t). (6.9)
∂t ∂x
Si además suponemos que no hay fuentes de calor (Q = 0), dividiendo la ecuación por cρ
la ecuación en derivadas parciales se vuelve
∂u ∂ 2u ∂u ∂ 2u
(x, t) = k 2 (x, t), ó = k 2, (6.10)
∂t ∂x ∂t ∂x
K0
donde la constante k = se llama difusividad térmica, la conductividad térmica divi-
cρ
dida por el producto del calor especı́fico y la densidad de masa. La ecuación (6.10) se
conoce como ecuación del calor ; corresponde a un caso sin fuentes de calor y propiedades
térmicas constantes. Si la energı́a calórica está inicialmente concentrada en un lugar, la
ecuación (6.10) describe cómo se va a propagar la energı́a, un proceso fı́sico conocido co-
mo difusión. Además de la temperatura, existen otras cantidades fı́sicas que se propagan
de forma similar satisfaciendo la misma ecuación diferencial (6.10). Por este motivo, la
ecuación (6.10) se conoce como ecuación de difusión. Por ejemplo, la concentración u(x, t)
de quı́micos (como ser perfumes y contaminantes) satisface la ecuación de difusión (6.10)
en ciertas situaciones unidimensionales.
u(x, 0) = f (x), 0 ≤ x ≤ L.
104
extremos x = 0 y x = L. Sin esta información, no podemos predecir el futuro. Como
tenemos dos derivadas con respecto a x en la ecuación diferencial, necesitaremos dos
condiciones en los extremos de la barra, generalmente una condición en cada extremo.
Discutiremos estas condiciones, conocidas como condiciones de borde (CB) en la siguiente
sección.
Frontera aislada. En otras situaciones es posible suponer que se conoce (o está pres-
cripto) el flujo de calor en lugar de la temperatura, es decir
∂u
− K0 (0, t) = φ(t), (6.12)
∂x
donde φ(t) es una función conocida o dada. Esto es equivalente a dar una condición para
la primera derivada ∂u/∂x en x = 0, es decir, se considera conocida o predeterminada la
pendiente de la curva y = u(x, t). La ecuación (6.12) no puede ser integrada en x puesto
que la pendiente se conoce sólo en el punto x = 0. El ejemplo más simple de flujo de
calor prescripto en la frontera es cuando uno de los extremos está perfectamente aislado.
En este caso no hay flujo en la frontera. Si el extremo x = 0 está aislado, luego
∂u
(0, t) = 0. (6.13)
∂x
105
muestran que el flujo de calor que sale de la barra es proporcional a la diferencia de tem-
peratura entre la barra y la temperatura del fluido exterior. Esta condición de frontera
se llama Ley de enfriamiento de Newton. Si la queremos escribir en el extremo x = 0, es
de la forma
∂u
− K0 (0) (0, t) = −H[u(0, t) − uB (t)], (6.14)
∂x
donde la constante de proporcionalidad H se lama coeficiente de transferencia de calor o
coeficiente de convección. Esta condición de borde involucra una combinación lineal de u y
∂u
. Debemos ser cuidadosos con el signo de la constante H de proporcionalidad. Veamos,
∂x
si la barra está más caliente que el baño del fluido (u(0, t) > uB (t)), entonces el calor
fluirá hacia fuera de la barra en x = 0. Luego, el calor fluye hacia la izquierda y en este
caso el flujo φ(x, t) debe ser negativo. Este es el motivo por el que aparece el signo menos
adelante de H en (6.14). Si hubiéramos supuesto u(0, t) < uB (t) habrı́amos arribado a
la misma conclusión respecto al signo frente a H. En un ejercicio se pide repetir este
razonamiento, para obtener que la ley de enfriamiento de Newton en el extremo derecho
de la barra (x = L) es
∂u
− K0 (L) (L, t) = H[u(L, t) − uB (t)], (6.15)
∂x
donde uB (t) es la temperatura externa en x = L. Inmediatamente puede apreciarse la di-
ferencia importante de signo entre la frontera izquierda (6.14) y la frontera derecha (6.15).
El coeficiente H en la ley de enfriamiento de Newton se determina experimentalmente.
Depende de las propiedades de la barra como ası́ también de las del fluido (inclusive de
la velocidad con la que el fluido se mueve). Si el coeficiente es muy pequeño, muy poca
energı́a fluye a través del extremo. En el lı́mite H → 0, la ley de enfriamiento de Newton,
se convierte en la condición de borde aislado. Podemos pensar que la ley de Newton con
H 6= 0 representa un extremo que no está perfectamente aislado.
Cuando H es grande, mucha energı́a fluye a través del extremo. En el lı́mite H → ∞,
la condición de borde se transforma en la condición de temperatura prescripta u(0, t) =
uB (t). Esto puede verse fácilmente si dividimos (6.14) por H:
K0 (0) ∂u
− (0, t) = −[u(0, t) − uB (t)],
H ∂x
y tomamos lı́mite cuando H → ∞.
106
Una condición de borde ocurre en cada extremo. No es necesario que las condiciones
en los dos extremos sean del mismo tipo. Por ejemplo, es posible que para x = 0 se tenga
una temperatura prescripta oscilante
y una condición de borde en cada extremo. Por ejemplo, cada extremo puede estar en
contacto con grandes baños de fluidos diferentes, tal que la temperatura en cada extremo
está prescripta
u(0, t) = T1 (t) u(L, t) = T2 (t). (6.18)
Uno de los objetivos de este curso es aprender a resolver el problema (6.16)–(6.18).
Tratemos ahora de ver si es posible hallar una solución estacionaria. Supongamos que
∂u
u es una solución en estado estacionario, entonces u no depende de t, y luego (x) = 0.
∂t
107
∂ 2u
Como u es solución de (6.16), se tiene que k = 0, pero como u = u(x) depende de
∂x2
una sola variable, las derivadas no son parciales sino ordinarias y tenemos que u = u(x)
satisface la ecuación diferencial ordinaria
d2 u
= 0. (6.20)
dx2
u(x) = C1 x + C2 .
Recordemos que T1 y T2 son datos dados y que tendrán un cierto valor en cada caso
particular, y que las incógnitas son C1 y C2 . Resolviendo (6.21) obtenemos que C2 = T1
(C2 = 10) y C1 = (T2 − T1 )/L (C1 = (12 − 10)/L = 2/L). Por lo tanto la única solución
en estado de equilibrio de la ecuación del calor con condiciones de temperatura prescripta
en los bordes es
T2 − T1 2
u(x) = T1 + x (u(x) = 10 + x).
L L
Si la barra fuera por ejemplo de longitud 4, entonces u(x) = 10 + x/2.
108
Convergencia al estado de equilibrio. Para el problema dependiente del tiem-
po (6.16) y (6.17) con condiciones de borde estacionarias (6.19), se espera que la dis-
tribución de temperatura se modifique a medida que pase el tiempo; no va a permanecer
igual a la distribución inicial f (x) a menos que ésta coincida con la solución estacionaria.
Si esperamos un perı́odo largo de tiempo, imaginamos que la influencia de las condiciones
en los extremos va a dominar el comportamiento de la temperatura en la barra: la con-
dición inicial poco importa si estamos mirando lo que ocurre después de mucho tiempo.
Se espera que eventualmente el perfil de la temperatura se aproxime a la distribución de
la temperatura en equilibrio
T2 − T1
lı́m u(x, t) = u(x) = T1 + x.
t→∞ L
Más adelante veremos cómo resolver el problema evolutivo (ecuación dependiente del
tiempo) y veremos efectivamente que este lı́mite es cierto. Mientras tanto, si sólo nos
interesa encontrar el estado estacionario o lı́mite cuando t → ∞, basta con resolver el
problema estacionario que hemos visto en el párrafo anterior.
Conclusión: Para hallar la solución estacionaria a la ecuación del calor, basta con
reemplazar ∂u/∂t por 0, y resolver la ecuación diferencial ordinaria en la variable x
resultante.
∂u ∂ 2u
EDP: =k 2 (6.22)
∂t ∂x
CI: u(x, 0) = f (x) (6.23)
∂u
CB1: (0, t) = 0 (6.24)
∂x
∂u
CB2: (L, t) = 0. (6.25)
∂x
∂u
El problema estacionario correspondiente se obtiene haciendo = 0. Luego la distribu-
∂t
ción de temperatura en estado de equilibrio satisface
d2 u
EDO: =0 (6.26)
dx2
du
CB1: (0) = 0 (6.27)
dx
du
CB2: (L) = 0, (6.28)
dx
109
d2 u
donde la condición inicial es ignorada (por el momento). La solución general de 2 = 0
dx
es otra vez una lı́nea recta arbitraria
u(x) = C1 x + C2 .
Las condiciones de borde implican que la pendiente debe ser cero en ambos extremos.
Geométricamente, cualquier lı́nea recta horizontal satisfará ambas condiciones de borde
simultáneamente. Luego la solución es cualquier distribución de temperatura constante,
du
por ejemplo: u(x) = 3, u(x) = 25.7. Algebraicamente, = C1 y entonces las dos
dx
condiciones de borde implican C1 = 0. Luego
u(x) = C2
d
E(t) = 0,
dt
110
por lo que E(t) es constante. En consecuencia la energı́a térmica total inicial E(0) debe
ser igual a la energı́a térmica total “final” E(∞) = lı́mt→∞ E(t). Pero
Z L Z L
E(0) = cρ u(0, t) dx = cρ f (x) dx
0 0
y
Z L Z L Z L
E(∞) = lı́m E(t) = lı́m cρ u(x, t) dx = cρ lı́m u(x, t) dx = cρ C2 dx.
t→∞ t→∞ 0 0 t→∞ 0
1 L
Z
u(x) = C2 = f (x) dx.
L 0
Esto es, la solución en estado de equilibrio para la ecuación del calor con bordes aislados
es la distribución de temperatura constante cuyo valor es el promedio de la distribución
de temperatura inicial. En este caso (bordes aislados) la condición inicial es importante, y
la barra no se olvida completamente de la condición inicial, “se acuerda de su promedio.”
En el caso en que no hay fuentes de calor Q = 0 y las propiedades térmicas son constantes,
la ecuación en derivadas parciales luce
∂u ∂ 2u
= k 2,
∂t ∂x
donde k = K0 /cρ. Antes de comenzar a estudiar métodos para resolver esta ecuación,
formularemos las ecuaciones en derivadas parciales que corresponden a problemas de flujo
de calor en dos o tres dimensiones espaciales. Veremos que la derivación y obtención de la
ecuación es muy similar a la que hemos hecho para el problema unidimensional, aunque
aparecerán diferencias importantes. Nos proponemos derivar ecuaciones nuevas y más
complejas (antes de resolver las más simples) para que cuando discutamos técnicas para
la resolución de las EDPs, tengamos más de un ejemplo para trabajar.
111
Energı́a calórica. Supongamos ahora que tenemos un sólido sobre el que queremos
estudiar la propagación de energı́a calórica. Comenzamos nuestra derivación de las ecua-
ciones o modelo matemático considerando una sub-región arbitraria R contenida en el
sólido. Como en el caso unidimensional, la conservación de la energı́a calórica se resume
en la siguiente ecuación:
razón de cambio energı́a calórica que
energı́a generada
de la energı́a fluye por la frontera
= + dentro de la región R
en la región R de la región R
por unidad de tiempo.
con respecto al tiempo por unidad de tiempo
Como ahora estamos considerando el caso de una región tridimensional, las constantes
que dependen del material dependerán de tres variables: c = c(x, y, z), ρ = ρ(x, y, z); y la
temperatura dependerá de cuatro variables, tres (x, y y z) que dan la posición espacial y
la cuarta (t) que representa el tiempo: u = u(x, y, z; t). Ahora bien, la energı́a dentro de
una sub-región arbitraria R es la integral triple
ZZZ
energı́a calórica en R = cρu dV.
R
De aquı́ en adelante, dV significa integral de volumen con respecto a las tres variables x,
y, z: dV = dx dy dz.
Flujo de calor y vectores normales. Ahora necesitamos una expresión para el flujo
de energı́a. En el problema unidimensional el flujo de calor φ se define como la cantidad
de calor por unidad de área que fluye hacia la derecha (φ < 0 significa que el flujo es
hacia la izquierda). En un problema tridimensional, el flujo de calor se dirige en alguna
dirección (no necesariamente alineada con alguno de los ejes), por lo que el flujo es un
campo vectorial φ. La magnitud de φ es la cantidad de energı́a calórica que fluye por
unidad de tiempo por unidad de área.
Por otro lado, cuando se considera la conservación de la energı́a, lo único que interesa
del flujo φ es la energı́a que pasa por unidad de tiempo a través de la superficie que
encierra la región R (denotaremos con ∂R a esta superficie que también se denomina
frontera de R). Es decir, lo que interesa de φ es la componente normal a la superficie
∂R. En cada punto de ∂R hay dos normales, la normal exterior, (que denotaremos con
n) y la normal interior (−n).
112
El signo de esta integral tiene la siguiente interpretación fı́sica: si es positiva, luego la
energı́a calórica total dentro de la barra decrece, y si es negativa, crece.
Denotemos con Q la densidad de energı́a térmica generada por unidad de tiempo por
unidad de volumen,
RRR luego el total de energı́a térmica generada por unidad de tiempo en
la región R es R
Q dV .
En consecuencia, la ley de conservación de energı́a calórica para una región arbitraria
tridimensional R se escribe
ZZZ ZZ ZZZ
d
cρu dV = − φ · n dσ + Q dV. (6.29)
dt R ∂R R
es decir, el flujo a través de los extremos puede expresarse como una integral sobre toda
la porción de la barra en consideración. El teorema de la divergencia nos permite realizar
algo análogo para tres variables. En efecto, el teorema de la divergencia dice que si F es
un campo vectorial de clase C 1 en una región que contiene a la región R, entonces
ZZZ ZZ
div F dV = F · n dσ.
R ∂R
Recordemos aquı́ que φ = φ(x, y, z; t) y estamos considerando siempre t fijo, por lo que
la divergencia se toma con respecto a las tres variables espaciales x, y, z. Es decir, si
φ(x, y, z; t) = φ1 (x, y, z; t)i + φ2 (x, y, z; t)j + φ3 (x, y, z; t)k,
113
Como la región R bajo consideración se encuentra fija, podemos pasar la derivada tem-
poral dentro del signo de la integral, y juntando todos los términos del mismo lado
obtenemos ZZZ
∂u
cρ + div φ − Q dV = 0.
R ∂t
Ahora bien, como esta integral es cero para toda sub-región R que uno quiera considerar,
es necesario que el integrando sea cero en todo punto:
∂u
cρ + div φ − Q = 0,
∂t
o equivalentemente
∂u
cρ = − div φ + Q. (6.31)
∂t
u(x, y, z; t) = T (x, y, z; t)
en todo punto de la frontera, donde T es una función conocida (dato) definida para todo
t en cada punto (x, y, z) de la frontera. También es posible que el flujo a través de la
frontera esté prescripto. Frecuentemente consideraremos el caso de frontera (o parte de
ella) aislada. Esto quiere decir que no hay flujo de calor a través de esa porción de la
frontera. Como el vector de flujo de calor es −K0 ∇u, el calor que fluye hacia fuera de
la superficie por unidad de área es la componente normal del vector de flujo, es decir
−K0 ∇u · n, donde n es el vector normal unitario exterior. Por lo tanto, en los puntos de
la frontera que se encuentran térmicamente aislados se tiene que
∇u · n = 0.
−K0 ∇u · n = H(u − uB ).
0 = div(K0 ∇u) + Q.
Notemos ahora que al estar en dimensión tres, la ecuación de estado estacionario sigue
siendo una ecuación en derivadas parciales. Es cierto que tenemos una variable menos (t
no aparece en esta ecuación), pero sigue habiendo derivadas con respecto a x, y y z. En
115
el caso que las propiedades térmicas sean constantes, la distribución de temperatura en
estado de equilibrio satisfará la ecuación
Q
∇2 u = − ,
K0
que se conoce como ecuación de Poisson.
Si además no hubiera fuentes de calor (Q = 0), la ecuación resultarı́a
∇2 u = 0;
6.6. Ejercicios
R1
6.1. Hallar una función continua en [0, 1] tal que 0
f (x) dx = 0 pero f no sea la función
constantemente cero.
6.2. Si se conoce u(x, t) ρ(x), c(x), K0 (x), dar una expresión para la energı́a térmica
total contenida en una barra de longitud L.
6.3. Considerar una barra unidimensional 0 ≤ x ≤ L. Suponer que la energı́a que fluye
hacia afuera de la barra en x = L es proporcional a la diferencia de temperatura entre
el extremo derecho de la barra y una temperatura externa conocida. Escribir la ley de
enfriamiento de Newton en términos matemáticos usando la temperatura u(x, t) y la
temperatura exterior uB (t). Analizar si los signos están correctos considerando el caso
u(L, t) > uB (t) y el caso u(L, t) < uB (t).
(b) No se pierde nada de energı́a calórica en x = x0 (es decir, la energı́a calórica que
fluye hacia afuera de una barra fluye hacia adentro de la otra).
116
Escribir ecuaciones matemáticas que representen estas dos condiciones. Para ello, denotar
con c1 , c2 , ρ1 , ρ2 , K01 , K02 las capacidades calóricas, las densidades y las conductividades
térmicas de ambos materiales.
6.5. Determinar la distribución de temperatura en estado de equilibrio para una barra
unidimensional con propiedades térmicas constantes con las siguientes fuentes y condi-
ciones de borde (suponer u(x, 0) = f (x))
(b) Determinar la energı́a calórica fluyendo hacia fuera de la barra por unidad de tiempo
en x = 0 y x = L.
(c) ¿Qué relación deberı́a existir entre las respuestas a la parte (a) y (b)?
6.7. Determinar la distribución de temperatura en equilibrio para una barra unidimen-
sional compuesta de dos materiales diferentes en contacto térmico perfecto en x = 1.
Para 0 ≤ x < 1, hay un material con coeficientes cρ = 1 y K0 = 1 con una fuente de
calor constante Q = 1. Para la otra parte de la barra 1 < x ≤ 2 no hay fuentes (Q = 0)
y los coeficientes son cρ = 2, K0 = 2, con condiciones de borde u(0) = 0 y u(2) = 0.
6.8. Si ambos extremos están aislados, concluir a partir de la ecuación en derivadas
parciales
∂u ∂ ∂u
c(x)ρ(x) (x, t) = (K0 (x) (x, t)), 0 ≤ x ≤ L, t ≥ 0,
∂t ∂x ∂x
que la energı́a térmica total en la barra es constante.
117
6.9. Considerar una barra unidimensional 0 ≤ x ≤ L de longitud conocida y propiedades
térmicas conocidas constantes y sin fuentes de calor, en estado de equilibrio. Suponer que
la temperatura es una constante desconocida T en x = L. Determinar T si se conocen la
temperatura y el flujo en x = 0. (Llamar T0 a la temperatura y φ0 al flujo en x = 0, y
encontrar una fórmula para T en términos de c, ρ, K0 , T0 y φ0 ).
6.10. Los dos extremos de una barra uniforme de longitud L están aislados. Hay una
fuente constante de energı́a térmica Q0 6= 0 y la temperatura es inicialmente u(x, 0) =
f (x).
(b) Calcular la energı́a térmica total a tiempo t. Dar una expresión en términos de f (x)
y de Q0 .
6.11. Para los siguientes problemas, determinar una distribución de temperatura en esta-
do de equilibrio (si existe). ¿Para qué valores de β existen soluciones? Explicar fı́sicamente.
∂u ∂ 2u ∂u ∂u
(a) = + 1, u(x, 0) = f (x) (0, t) = 1 (L, t) = β
∂t ∂x2 ∂x ∂x
∂u ∂ 2u ∂u ∂u
(b) = , u(x, 0) = f (x) (0, t) = 1 (L, t) = β
∂t ∂x2 ∂x ∂x
∂u ∂ 2u ∂u ∂u
(c) = + x − β, u(x, 0) = f (x) (0, t) = 0 (L, t) = 0
∂t ∂x2 ∂x ∂x
6.12. Denotemos con c(x, y, z; t) la concentración de un contaminante (cantidad de con-
taminante por unidad de volumen)
(c) Derivar la ecuación en derivadas parciales que gobierna la difusión del contaminante.
para toda superficie cerrada S. (Ayuda: usar el teorema de la divergencia). Dar una
interpretación fı́sica de este resultado (en el contexto de flujo de calor).
118
Capı́tulo 7
y dos condiciones de borde. Por ejemplo, si ambos extremos de la barra tienen la tempe-
ratura prescripta, entonces
7.1. Linealidad
Ası́ como en el estudio de ecuaciones diferenciales ordinarias, el concepto de linealidad
será para nosotros de gran importancia. Un operador L es lineal cuando para todo par
de funciones u1 y u2 , y todo par de números reales c1 , c2 se satisface
119
Otros ejemplos de operadores lineales son ∂/∂t, y ∂ 2 /∂x2 . Y otro un poquito más com-
plicado es el operador conocido como operador del calor:
∂ ∂2
− k 2,
∂t ∂x
que también es un operador lineal.
Una ecuación se dice lineal si es de la forma
L(u) = f, (7.5)
con L un operador lineal y f conocida (dato). Los siguientes son ejemplos de ecuaciones
en derivadas parciales lineales
∂u ∂ 2u
=k + f (x, t) (7.6)
∂t ∂x2
∂u ∂ 2u
= k 2 + α(x, t)u + f (x, t) (7.7)
∂t ∂x
∂ 2u ∂ 2u
+ =0 (7.8)
∂x2 ∂y 2
∂u ∂ 3u
= + α(x, t)u. (7.9)
∂t ∂x3
Los siguientes son ejemplos de ecuaciones en derivadas parciales no-lineales
∂u ∂ 2u
= k 2 + α(x, t)u4 (7.10)
∂t ∂x
∂u ∂u ∂ 3u
+u = . (7.11)
∂t ∂x ∂x3
∂u
Los términos u4 y u son no-lineales; no satisfacen (7.4).
∂x
Si f = 0, entonces (7.5) se llama ecuación lineal homogénea. La ecuación del calor
∂u ∂ 2u
−k 2 =0 (7.12)
∂t ∂x
es un ejemplo de ecuación lineal homogénea, como también lo son (7.8) y (7.9). De (7.4),
tomando c1 = c2 = 0, se sigue que si aplicamos un operador lineal L a la función u ≡ 0,
obtenemos L(0) = 0. Por lo tanto, u ≡ 0 es siempre solución de una ecuación lineal
homogénea. Por ejemplo, u ≡ 0 satisface la ecuación del calor (7.12). La solución u ≡ 0 se
llama solución trivial. La manera más simple de ver si una ecuación es homogénea consiste
en reemplazar u por la función idénticamente cero. Si u ≡ 0 satisface una ecuación lineal,
entonces f = 0 y la ecuación lineal es homogénea.
La propiedad fundamental de los operadores lineales es que nos permiten sumar las
soluciones de una misma ecuación lineal y obtener otra solución:
120
La demostración de este principio depende de la definición de operador lineal. Supon-
gamos que u1 y u2 son dos soluciones de la ecuación lineal homogénea L(u) = 0, es decir,
L(u1 ) = 0 y L(u2 ) = 0. Calculemos L(u) para u = c1 u1 + c2 u2 . De la definición (7.4) de
operador lineal
∂u ∂ 2u
EDP: =k 2 0≤x≤L t>0 (7.18)
∂t
( ∂x
u(0, t) = 0
CB: t≥0 (7.19)
u(L, t) = 0
CI: u(x, 0) = f (x) 0 ≤ x ≤ L. (7.20)
El problema consiste en una ecuación en derivadas parciales lineal homogénea con con-
diciones de borde lineales homogéneas. Hay dos razones para investigar este tipo de
problemas, más allá de que lo podremos resolver por el método de separación de varia-
bles. Primero, este problema es relevante desde el punto de vista fı́sico, correspondiente
a una barra unidimensional sin fuentes y con ambos extremos inmersos en un baño a
temperatura 0o . Estamos interesado en poder predecir cómo cambia la energı́a térmica
121
inicial (representada en el dato inicial f (x)) a medida que pasa el tiempo en esta situación
fı́sica relativamente sencilla. Segundo, veremos que para poder resolver el problema no-
homogéneo (7.1)–(7.3), necesitaremos saber cómo resolver el problema homogéneo (7.18)–
(7.20).
Buscamos todas las soluciones posibles que sean producto de una función de x y una
de t y que satisfagan la ecuación diferencial y las condiciones de borde homogéneas.
∂u dG(t) ∂ 2u d2 φ(x)
(x, t) = φ(x) , (x, t) = G(t).
∂t dt ∂x2 dx2
En consecuencia la ecuación del calor (7.18) implica que
dG(t) d2 φ(x)
φ(x) =k G(t). (7.22)
dt dx2
Ahora nos damos cuenta de que podemos “separar variables” dividiendo ambos miembros
de (7.22) por el producto φ(x)G(t):
1 dG(t) 1 d2 φ(x)
=k .
G(t) dt φ(x) dx2
Se ve que las variables han sido “separadas” en el sentido que el lado izquierdo es sólo
función de la variable t y el lado derecho es sólo función de la variable x. Podemos
122
continuar trabajando con esta expresión, pero es conveniente dividir también por k, y
luego
1 dG 1 d2 φ
= . (7.23)
kG dt
| {z } φ dx2
| {z }
función de t función de x
¿Cómo es posible que una función solamente dependiente de t sea igual a una función que
depende solamente de x? La única posibilidad es que ambas expresiones sean constantes
e independientes tanto de x como de t. Y por la igualdad deben ser iguales a la misma
constante:
1 dG 1 d2 φ
= = −λ, (7.24)
kG dt φ dx2
d2 φ
= −λφ (7.25)
dx2
dG
= −λkG. (7.26)
dt
Tenemos ahora tres incógnitas: φ, G, y λ. Para hallar las soluciones producto (7.21)
debemos resolver las dos equaciones diferenciales ordinarias que acabamos de deducir, y
ver para qué constantes λ estas soluciones no dan la solucion trivial.
Repetimos que λ es constante, y es la misma en ambas ecuaciones (7.25) y (7.24). Las
soluciones producto también deben satisfacer las dos condiciones de borde. Por ejemplo,
u(0, t) = 0 implica que φ(0)G(t) = 0. Hay aquı́ dos posibilidades. O bien G(t) ≡ 0 (≡
quiere decir que es idénticamente cero, para todo t), o φ(0) = 0. Si G(t) ≡ 0, entonces
u ≡ 0 y obtuvimos la solución trivial. Como estamos buscando soluciones no-triviales,
descartamos este caso y suponemos G(t) 6≡ 0. Si G(t) 6≡ 0, es necesario que
φ(0) = 0. (7.27)
φ(L) = 0. (7.28)
123
Conclusión: para encontrar las soluciones de variables separables u(x, t) =
φ(x)G(t), debemos encontrar las ternas (φ, G, λ) para las cuales se satisfagan las
siguientes ecuaciones simultáneamente:
d2 φ
= −λφ
dx2
dG
= −λkG.
dt
d2 φ
= −λφ
dx2
(7.30)
φ(0) = 0
φ(L) = 0.
124
No hay ninguna teorı́a simple que garantice que la solución existe o es única para este tipo
de problemas. En particular, notamos que φ(x) ≡ 0 satisface la EDO y ambas condiciones
de borde homogéneas, sin importar el valor de λ, aún para λ < 0. Esta solución se conoce
con el nombre de solución trivial del problema a valores de borde. Como u(x, t) =
φ(x)G(t), corresponde a u(x, t) ≡ 0, si las soluciones de (7.30) fueran únicas, φ(x) ≡ 0
serı́a la única solución; y no podrı́amos obtener soluciones no-triviales de una EDP lineal
homogénea por el método de separación de variables. Afortunadamente sı́ existen otras
soluciones de (7.30), aunque no para cualquier valor de λ sino sólo para algunos. Veremos
a continuación que hay ciertos valores de λ, llamados autovalores de (7.30), para los
cuales existen soluciones no-triviales φ(x). Una φ(x) no-trivial, correspondiente a algún
valor de λ se llama autofunción correspondiente al autovalor λ.
1. λ >
√ 0: las dos raı́ces son puramente imaginarias y son los números complejos r =
±i λ.
125
Si c2 = 0, obtenemos la solución trivial φ(x) ≡ 0 que no nos interesa. Si c2 6= 0 entonces
tenemos que buscar todos los valores de λ para los que
√
sen λL = 0.
√
En√ otras palabras, λL debe ser un cero de la función seno, es decir un múltiplo entero de
π: λL = nπ. Como hemos supuesto λ > 0, n toma los valores 1, 2, . . . . Los autovalores
λ son entonces nπ 2
λ= , n = 1, 2, 3 . . . . (7.33)
L
2
La autofunción correspondiente al autovalor λ = nπ L
es
√ nπx
φ(x) = c2 sen λx = c2 sen (7.34)
L
donde c2 es una constante arbitraria. A menudo elegimos una autofunción seleccionando
un valor no-nulo para c2 , por ejemplo c2 = 1. Debemos recordar que cualquier autofunción
puede ser multiplicada por una constante arbitraria y seguir siendo autofunción del
mismo autovalor pues la ecuación es lineal y homogénea.
127
Por lo tanto, todas las soluciones producto de la ecuación del calor son
nπx −( nπ 2
u(x, t) = B sen e L ) kt , n = 1, 2, 3, . . .
L
donde B es una constante arbitraria. Es decir, para cada número real B y cada entero
positivo n, la expresión anterior es una solución de la ecuación del calor con condiciones
de borde de temperatura nula.
Ejemplos: las siguientes funciones son soluciones de la ecuación del calor con condi-
ciones de borde de temperatura nula:
3πx − 9π22 kt
u1 (x, t) = sen e L
L
πx − π22 kt
u2 (x, t) = 4.5 sen e L .
L
Y como la ecuación y las condiciones de borde que estamos considerando son lineales y
homogéneas, también
3πx − 9π22 kt πx − π22 kt
u3 (x, t) = sen e L + 4.5 sen e L
L L
es una solución.
¿Qué condición inicial satisfacen estas tres soluciones? Basta reemplazar t por cero
para ver que
3πx
u1 (x, 0) = sen
L
πx
u2 (x, 0) = 4.5 sen
L
3πx πx
u3 (x, 0) = sen + 4.5 sen .
L L
No parece haber mucha libertad para elegir condiciones iniciales generales. Nos ocupare-
mos de este tema en la próxima sección.
128
Principio de superposición. Las soluciones producto son soluciones muy particu-
lares, pues sólo sirven cuando la condición inicial tiene la forma apropiada u(x, 0) =
B sen nπx
L
. Sin embargo, este tipo de soluciones es muy útil también en muchas otras
situaciones. Consideremos la misma EDP con las mismas condiciones de borde pero con
la condición inicial
3πx 8πx
u(x, 0) = 4 sen + 7 sen .
L L
La solución de este problema puede obtenerse sumando dos soluciones producto:
3πx −( 3π 2 8πx −( 8π 2
u(x, t) = 4 sen e L ) kt + 7 sen e L ) kt .
L L
Inmediatamente vemos esta función resuelve la EDP, las CB y también la CI.
es también una solución. Como sabemos del método de separación de variables que la
nπ 2
función sen nπx e−( L ) kt es solución de la ecuación del calor para todo entero n positivo,
L
entonces cualquier combinación lineal de estas soluciones es también solución de la misma
ecuación del calor con las mismas condiciones de borde homogéneas. Por lo tanto
M
nπx −( nπ 2
e L ) kt
X
u(x, t) = Bn sen (7.37)
n=1
L
es solución de la ecuación del calor (con condiciones de borde nulas) para cualquier M y
cualesquiera valores de Bn . Entonces hemos ahora extendido la familia de soluciones. De
soluciones producto, hemos pasado a combinaciones lineales de soluciones produc-
to. Si en (7.37) tomamos t = 0 vemos que
M
X nπx
u(x, 0) = Bn sen .
n=1
L
Por lo que podemos resolver la ecuación del calor para cualquier condición inicial de la
forma
M
X nπx
u(x, 0) = f (x) = Bn sen .
n=1
L
¿Qué hacemos si tenemos que resolver la ecuación del calor con un dato inicial f (x) que
no es de esta forma? La teorı́a de las series de Fourier dice que:
1. Cualquier función f (x) (con algunas restricciones razonables) puede ser aproximada
(en algún sentido) por combinaciones lineales de sen nπx
L
.
129
2. La aproximación tal vez no será buena para M pequeño, pero mejora a medida que
M crece.
Afirmamos entonces que cualquier dato inicial f (x) se puede escribir como una com-
binación lineal infinita de sen nπx
L
, conocida como serie de Fourier:
∞
X nπx
f (x) = Bn sen . (7.38)
n=1
L
Volvamos a lo que nos interesa, utilizar esta propiedad para encontrar un modo de deter-
minar Bn . Supongamos entonces que
∞
X nπx
f (x) = Bn sen ,
n=1
L
130
multipliquemos ambos lados por sen mπx
L
para un m fijo independiente de n. Obtenemos
∞
mπx X nπx mπx
f (x) sen = Bn sen sen .
L n=1
L L
Ahora bien, ¿a qué equivale cada término de la suma infinita del lado derecho? Recor-
demos que m está fijo (por ejemplo, supongamos m = 23), y que queremos evaluar las
integrales para cada n desde 1 hasta ∞. Entonces el primer término es
Z L (
1πx mπx 0 si m 6= 1
sen sen dx =
0 L L L/2 si m = 1
( = 0 si m = 23.)
7.2.7. Resumen
Resumamos el método de separación de variables para el ejemplo
∂u ∂ 2u
EDP: =k 2 t ≥ 0, 0≤x≤L
∂t ∂x
CB1: u(0, t) = 0 t≥0
CB2: u(L, t) = 0 t≥0
CI: u(x, 0) = f (x) 0 ≤ x ≤ L.
1. Asegurarse de tener una EDP lineal y homogénea con condiciones de borde homogéneas.
131
4. Determinar las constantes de separación como los autovalores de un problema a valores
de borde (utilizar las condiciones de borde CB).
5. Resolver las otras ecuaciones diferenciales ordinarias para los autovalores encontrados.
Indicar todas las soluciones producto de la EDP que se obtienen por este método.
6. Aplicar el principio de superposición (para una combinación lineal de todas las solu-
ciones producto)
8. Determinar los coeficientes utilizando alguna relación integral entre las autofunciones.
∂u ∂ 2u
EDP: =k 2 t ≥ 0, 0≤x≤L
∂t ∂x
∂u
CB1: (0, t) = 0 t≥0
∂x
∂u
CB2: (L, t) = 0 t≥0
∂x
CI: u(x, 0) = f (x) 0 ≤ x ≤ L.
Como repaso, mencionamos que es un problema de conducción del calor en una barra
unidimensional con propiedades térmicas constantes y sin fuentes. El problema es muy
similar al tratado en la Sección 7.2, la única diferencia reside en las condiciones de borde.
Aquı́ los bordes se encuentran aislados, mientras que en la Sección 7.2 los bordes estaban
mantenidos a temperatura 0o . Tanto la ecuación en derivadas parciales como las condi-
ciones de borde son lineales y homogéneas, por lo que podemos utilizar el método de
separación de variables.
Busquemos entonces todas las soluciones producto posibles, es decir, proponemos
u(x, t) = φ(x)G(t).
132
Si ésta fuera una solución, entonces deberı́a cumplirse que
dG d2 φ
= −λkG = −λφ,
dt dx2
donde λ es la constante de separación. Las condiciones de bordes aislados implican que
dφ dφ
(0)G(t) = 0 y (L)G(t) = 0 para todo t ≥ 0.
dx dx
dφ dφ
Para que la solución no sea la trivial, será necesario que (0) = 0 y (L) = 0. La
dx dx
constante de separación λ se determina entonces encontrando todos los λ para los que el
problema a valores en el borde:
d2 φ dφ dφ
= −λφ, (0) = 0, (L) = 0,
dx2 dx dx
tiene solución no-trivial. Aunque la ecuación diferencial ordinaria es la misma que la del
caso ya analizado, las condiciones de borde son diferentes. Dada esta diferencia, debemos
repetir el cálculo de los autovalores y autofunciones. Nuevamente debemos considerar tres
casos: λ > 0, λ = 0, λ < 0.
Analicemos primero el caso λ > 0: la solución general es
√ √
φ(x) = c1 cos λx + c2 sen λx.
Para que φ satisfaga las condiciones de frontera, debemos calcular
dφ √ √ √
(x) = λ −c1 sen λx + c2 cos λx .
dx
dφ √
La condición de borde (0) = 0 implica c2 λ = 0, y luego c2 = 0 pues λ > 0. Por lo
dx
tanto √
φ(x) = c1 cos λx.
Los autovalores λ y sus autofunciones correspondientes se determinan a partir de la
dφ
condición de borde restante (L) = 0:
dx
√ √
−c1 λ sen λL = 0.
√ √
Como antes, para soluciones no triviales c1 6= 0 y entonces sen λL = 0. Es decir λL =
nπ, ó bien nπ 2
λ= , n = 1, 2, 3, . . . .
L
Las autofunciones correspondientes son esta vez cosenos
nπx
φ(x) = c1 cos , n = 1, 2, 3 . . . .
L
Ahora bien, para cualquier valor de λ, las soluciones G(t) de la ecuación diferencial
dG
ordinaria = −λkG son las funciones de la forma
dt
G(t) = ce−λkt .
133
En consecuencia, las soluciones producto (para λ > 0) de la ecuación en derivadas par-
ciales son
nπx −( nπ 2
u(x, t) = A cos e L ) kt , (7.42)
L
para cualquier constante A, y para cualquier entero positivo n.
φ(x) = c1 + c2 x.
dφ
La derivada de φ es la función constante (x) = c2 , por lo que las dos condiciones de
dx
borde dan la misma condición c2 = 0. Por lo tanto, sı́ existen soluciones no triviales para
λ = 0. La autofunción corresondiente al autovalor λ = 0 es cualquier función constante
φ(x) = c1 .
La parte dependiente del tiempo G(t) también es constante, pues para λ = 0, la función
e−λkt ≡ e0 ≡ 1. Por lo tanto, otra solución producto de la ecuación del calor con bordes
aislados es
u(x, t) = A
para cualquier constante A.
Podemos analizar el caso λ < 0 desde un punto de vista fı́sico. Si hubiera algún
autovalor negativo, esto significarı́a que habrı́a alguna solución (no trivial) de la ecuación
de la forma u(x, t) = φ(x)e−λkt . Si λ < 0 entonces −λk > 0 y esta solución crece
exponencialmente a medida que pasa el tiempo. Pero esto es fı́sicamente imposible si los
bordes están aislados y no hay fuentes de calor (ambas hipótesis se cumplen en nuestro
caso). La contradicción provino de suponer la existencia de autovalores λ < 0. Concluimos
entonces que no hay autovalores negativos. En el ejercicio 7.3.4 se pide la demostración
matemática de que no hay autovalores negativos.
Las recién mencionadas son las soluciones producto de la ecuación diferencial que tam-
bién satisfacen las condiciones de borde, pero no satisfacen necesariamente la condición
inicial. Para satisfacer la condición inicial utilizamos el principio de superposición. Toma-
mos como solución general una combinación lineal de todas las soluciones producto que
hemos encontrado. Es decir
∞
nπx −( nπ 2
e L ) kt .
X
u(x, t) = A0 + An cos (7.43)
n=1
L
134
Es interesante notar que esta solución también se puede escribı́r ası́:
∞
nπx −( nπ 2
e L ) kt ,
X
u(x, t) = An cos
n=0
L
1 L
Z
A0 = f (x) dx
L 0
2 L
Z
mπx
Am = f (x) cos dx, m = 1, 2, 3 . . . .
L 0 L
Las dos fórmulas diferentes para los coeficientes Am son un poco molestas en esta serie de
cosenos. Sin embargo reflejan una diferencia significativa entre las soluciones producto de
135
la EDP para λ > 0 y las soluciones para λ = 0. Todas las soluciones para λ > 0 decaen
exponencialmente en el tiempo, mientras que la solución para λ = 0 permanece constante
en el tiempo. Por lo tanto, cuando t → ∞ la serie infinita (7.43) se aproxima al estado
estacionario
1 L
Z
lı́m u(x, t) = A0 = f (x) dx.
t→∞ L 0
No sólo la distribución de temperatura estacionaria es constante, A0 , sino que recono-
cemos la constante como el promedio de la distribución inicial de temperatura. Esto
coincide con las conclusiones que habı́amos obtenido al estudiar soluciones estacionarias
en la Sección 6.4.
Supongamos que un alambre delgado se curva para formar una circunferencia, como
en la Figura 7.1. Por razones que no son claras ahora, pero que lo serán más adelante,
supongamos que el alambre tiene longitud 2L, es decir, L es la mitad de la longitud del
2L L
alambre. Como la circunferencia tiene una longitud de 2πr, el radio es r = = .
2π π
Si el alambre es suficientemente delgado, es razonable suponer que la temperatura es
constante en cada sección transversal del alambre curvo. En esta situación la temperatura
del alambre debe satisfacer una ecuación del calor unidimensional, donde x representa la
longitud de arco a lo largo del alambre:
∂u ∂ 2u
= k 2. (7.46)
∂t ∂x
Nuevamente suponemos que el alambre tiene propiedades térmicas constantes y no hay
fuentes de calor. Es conveniente para simplificar algunos cálculos, suponer que la longitud
de arco x varı́a de −L a L (en lugar de variar de 0 a 2L).
Supongamos que el alambre está perfectamente conectado a sı́ mismo en los extremos
(x = −L a x = L). Allı́ deberı́an cumplirse las condiciones de contacto térmico perfecto:
La temperatura u(x, t) deberı́a ser continua
u(−L, t) = u(L, t). (7.47)
136
También el flujo de calor deberı́a ser continuo, y como la conductividad térmica es cons-
tante, esto implica que la derivada con respecto a x de la temperatura es continua:
∂u ∂u
(−L, t) = (L, t). (7.48)
∂x ∂x
Las dos condiciones de borde para la ecuación en derivadas parciales son entonces (7.47)
y (7.48). La condición inicial es una función dada
El problema matemático consiste luego de la EDP lineal homogénea (7.46), con condicio-
nes de borde lineales y homogéneas (7.47) y (7.48). Es por lo tanto plausible la utilización
del método de separación de variables.
Proponemos una solución producto u(x, t) = φ(x)G(t). Luego obtenemos para φ y G
las ecuaciones diferenciales ordinarias
dG d2 φ
= −λkG = −λφ.
dt dx2
Las condiciones de borde que nos interesan ahora implican que la función φ debe satisfacer
dφ dφ
φ(−L) = φ(L) y (−L) = (L). La constante de separación λ se determina entonces
dx dx
encontrando todos los λ para los que el problema a valores en el borde:
d2 φ dφ dφ
= −λφ, φ(−L) = φ(L) (−L) = (L),
dx2 dx dx
tiene solución no-trivial. Las condiciones de borde con que nos encontramos se denominan
comúnmente condiciones de borde periódicas porque aunque el problema se piensa
definido en −L < x < L, se puede pensar que está definido periódicamente para todo
x ∈ R; la temperatura será periódica (x = x0 es el mismo punto fı́sico que x = x0 + 2L,
y por lo tanto tendrá la misma temperatura). Si λ > 0, la solución general del problema
para φ es √ √
φ(x) = c1 cos λx + c2 sen λx.
La condición de borde φ(−L) = φ(L) implica que
√ √ √ √
c1 cos λ(−L) + c2 sen λ(−L) = c1 cos λL + c2 sen λL
√ √
Como
√ el coseno es una
√ función par, cos λ(−L) = cos λL, y como el seno es impar
sen λ(−L) = − sen λL. Se sigue entonces que φ(−L) = φ(L) se satisface cuando
√ √ √ √
c1 cos λL − c2 sen λL = c1 cos λL + c2 sen λL
Antes de resolver esta ecuación, analicemos la segunda condición de borde, que involucra
a la derivada,
dφ √ √ √
(x) = λ −c1 sen λx + c2 cos λx .
dx
137
dφ dφ
Luego (−L) = (L) se cumple si
dx dx
√ √ √ √ √ √
λ −c1 sen λ(−L) + c2 cos λ(−L) = λ −c1 sen λL + c2 cos λL
138
Conclusión: Las soluciones producto para la conducción del calor en un
alambre circular de longitud 2L son:
nπx −( nπ 2
uan (x, t) = an cos e L ) kt n = 1, 2, . . . ,
L
nπx −( nπ 2
b
un (x, t) = bn sen e L ) kt n = 1, 2, . . . ,
L
u0 (x, t) = a0 ,
Aquı́ la función f (x) es una combinación lineal de senos y de cosenos. Otra diferencia
crucial es que (7.53) es válida para −L ≤ x ≤ L, mientras que para las series de senos o
de cosenos de los ejemplos anteriores era válida para 0 ≤ x ≤ L.
Ahora queremos determinar los coeficientes a0 , an , bn (n ≥ 1) a partir de (7.53).
Nuevamente, las autofunciones satisfacen las siguientes propiedades integrales (notar que
ahora son diferentes pues estamos considerando el intervalo [−L, L] que es simétrico
respecto al origen):
Z L 0
si n 6= m
nπx mπx
cos cos dx = L si n = m 6= 0
−L L L
2L si n = m = 0
(
Z L
nπx mπx 0 si n 6= m (7.54)
sen sen dx =
−L L L L si n = m 6= 0
Z L
nπx mπx
y sen cos dx = 0 (siempre),
−L L L
139
Los coeficientes se obtienen de forma similar a como se hizo antes. Si multiplica-
mos (7.53) por cos mπx
L
y/o por sen mπx
L
(hacemos las dos cuentas a la vez) y luego inte-
gramos de x = −L a x = L, obtenemos
mπx mπx
cos L cos
Z L ∞ Z L
L
X nπx
f (x) dx = an cos dx
−L mπx n=0 −L L mπx
sen
sen
L L
mπx
∞ Z L cos
L
X nπx
+ bn sen . dx.
n=1 −L L mπx
sen
L
Ahora utilizamos las propiedades integrales (7.54) y obtenemos
Z L Z L
mπx mπx 2
f (x) cos dx = am cos dx
−L L −L L
Z L Z L
mπx mπx 2
f (x) sen dx = bm sen dx
−L L −L L
Despejando obtenemos
Z L
1
a0 = f (x) dx
2L −L
1 L
Z
mπx
am = f (x) cos dx
L −L L
1 L
Z
mπx
bm = f (x) sen dx
L −L L
donde las dos últimas fórmulas valen para m = 1, 2, 3, . . . .
140
Consideraremos primero el caso de condiciones de borde de temperatura prescripta:
CB1: u(0, y) = g1 (y) 0≤y≤H
CB2: u(L, y) = g2 (y) 0≤y≤H
CB3: u(x, 0) = f1 (x) 0≤x≤L
CB4: u(x, H) = f2 (x) 0≤x≤L
donde g1 (y), g2 (y), f1 (x) y f2 (x) son funciones dadas. Otras condiciones de borde serán
presentadas en los ejercicios. En este caso, la EDP es lineal y homogénea pero las condi-
ciones de borde, aunque lineales, no son homogéneas. No podremos aplicar el método de
separación de variables a este problema tal como está. La razón de esto es que al separar
variables, el problema a valores de borde que determina las constantes de separación po-
sibles o autovalores debe tener condiciones de borde homogéneas. En este ejemplo, todas
las condiciones de borde son no-homogéneas.
Podemos superar esta dificultad dándonos cuenta de que el problema original es no-
homogéneo debido a las cuatro condiciones de borde no-homogéneas. Utilizaremos el
principio de superposición para desglosar nuestro problema en cuatro subproblemas, cada
uno de los cuales tiene sólo una condición no-homogénea. Escribimos
u(x, t) = u1 (x, t) + u2 (x, t) + u3 (x, t) + u4 (x, t),
donde cada una de las funciones ui (x, y) (i = 1, 2, 3, 4) satisface la ecuación de Laplace con
una condición de borde no-homogénea y las otras tres condiciones de borde homogéneas,
como se muestra en la Figura 7.2.
El método para hallar cualquiera de las ui (x, y) es el mismo; sólo algunos detalles
difieren. Haremos la resolución para encontrar u4 (x, y), y dejaremos los otros casos para
141
los ejercicios:
∂ 2 u4 ∂ 2 u4
EDP: + =0 0 ≤ x ≤ L, 0≤y≤H
∂x2 ∂y 2
CB1: u4 (0, y) = g1 (y) 0≤y≤H
CB2: u4 (L, y) = 0 0≤y≤H
CB3: u4 (x, 0) = 0 0≤x≤L
CB4: u4 (x, H) = 0 0 ≤ x ≤ L.
Nos proponemos resolver este problema por el método de separación de variables. Comen-
zamos ignorando la condición no-homogénea u4 (0, y) = g1 (y). Eventualmente sumaremos
soluciones producto para obtener esta condición, pero por ahora la ignoraremos. Busca-
mos soluciones producto
u4 (x, y) = h(x)φ(y),
que satisfagan las condiciones de borde homogéneas CB2, CB3 y CB4. Éstas nos dicen
que
h(L) = 0, φ(0) = 0, φ(H) = 0.
Luego, la parte de la solución dependiente de y, φ(y) tiene dos condiciones de borde
homogéneas, mientras que la parte dependiente de x sólo tiene una. Si sustituimos la
solución producto en la ecuación de Laplace obtenemos
1 00 1 00
h (x) = − φ (y). (7.55)
h(x) φ(y)
El lado izquierdo es sólo función de la variable x, mientras que el lado derecho depende
sólo de y. Por lo tanto, ambos deben ser iguales a una constante de separación. ¿Qué
usamos? ¿λ ó −λ? Una de ellas es más conveniente. Si la constante de separación fuera
negativa (como lo era antes), la ecuación (7.55) implicarı́a que h(x) oscila (trigonométri-
ca) y φ(y) es una combinación de exponenciales. Pero las condiciones de borde para φ(y)
indican que este caso no es posible. Por otro lado, si la constante de separación fuera po-
sitiva, (7.55) implicarı́a que h(x) es combinación de exponenciales y φ(y) es combinación
de senos y cosenos. Esto parece más razonable dadas las condiciones de borde y entonces
introducimos la variable de separación λ (aunque no suponemos λ ≥ 0):
1 00 1 00
h (x) = − φ (y) = λ.
h(x) φ(y)
142
El problema dependiente de la variable x no es un problema a valores en la frontera pues
no tiene dos condiciones homogéneas:
nπ(x − L) nπ(x − L)
cosh y senh
H H
son soluciones linealmente independientes de (7.56). Por lo que la solución general también
puede escribirse como
nπ(x − L) nπ(x − L)
h(x) = a1 cosh + a2 senh .
H H
Resulta natural preguntarse ahora por qué es más conveniente esta expresión de la solu-
ción general. Bien, si queremos forzar la condición de borde h(L) = 0 obtenemos
nπ(L − L) nπ(L − L)
0 = h(L) = a1 cosh + a2 senh = a1 cosh 0 + a2 senh 0 = a1 .
H H
nπy nπ(x − L)
u4 (x, y) = A sen senh , n = 1, 2, . . . .
H H
143
Ahora queremos encontrar una solución más general, para lo cual utilizamos el prin-
cipio de superposición, obteniendo
∞
X nπy nπ(x − L)
u4 (x, y) = An sen senh .
n=1
H H
Una vez que tenemos esta solución general, imponemos la condición de borde no-homogénea
u4 (0, y) = g1 (y):
∞
X nπy nπ(−L)
u4 (0, y) = An sen senh = g1 (y).
n=1
H H
Este es el mismo tipo de serie de senos que ya hemos discutido con una sola diferencia:
Antes tenı́amos los coeficientes An multiplicando a las funciones seno, pero ahora tenemos
que An senh nπ(−L)
H
multiplica a las funciones seno (notar que An senh nπ(−L)
H
no depende
nπ(−L)
de la variable y). Por ende, si llamamos Bn = An senh H nuestro problema consiste
en hallar Bn tal que
∞
X nπy
Bn sen = g1 (y).
n=1
H
Luego, por lo que ya hemos visto en secciones anteriores
2 H
Z
nπy
Bn = g1 (y) sen dy,
H 0 H
y despejando An obtenemos
Z H
Bn 2 nπy
An = nπ(−L)
= nπ(−L)
g1 (y) sen dy.
senh H senh 0 H
H H
144
Figura 7.3: Ecuación de Lapla-
ce dentro de un disco
Luego, utilizando la regla de la cadena y la regla para la derivada del producto vemos
que
∂u ∂u ∂r ∂u ∂θ
= +
∂x ∂r ∂x ∂θ ∂x
∂ 2u
∂ ∂u ∂ ∂u ∂r ∂u ∂θ
= = +
∂x2 ∂x ∂x ∂x ∂r ∂x ∂θ ∂x
2 2
∂ 2 u ∂r ∂ 2 u ∂θ ∂r ∂u ∂ 2 r ∂ 2 u ∂r ∂θ ∂ 2 u ∂θ ∂u ∂ 2 θ
= 2 + + + + +
∂r ∂x ∂θ∂r ∂x ∂x ∂r ∂x2 ∂r∂θ ∂x ∂x ∂θ2 ∂x ∂θ ∂x2
2 2
∂ 2 u ∂r ∂u ∂ 2 r ∂ 2 u ∂r ∂θ ∂ 2 u ∂θ ∂u ∂ 2 θ
= 2 + + 2 + +
∂r ∂x ∂r ∂x2 ∂r∂θ ∂x ∂x ∂θ2 ∂x ∂θ ∂x2
Análogamente
2 2
∂ 2u ∂ 2u ∂u ∂ 2 r ∂ 2 u ∂r ∂θ ∂ 2 u ∂u ∂ 2 θ
∂r ∂θ
= + + 2 + 2 +
∂y 2 ∂r2 ∂y ∂r ∂y 2 ∂r∂θ ∂y ∂y ∂θ ∂y ∂θ ∂y 2
Utilizando ahora las expresiones para r(x, y) y θ(x, y), derivándolas y haciendo algo de
aritmética obtenemos (hacer estos pasos como ejercicio)
∂ 2u ∂ 2u ∂ 2 u 1 ∂u 1 ∂ 2u
∇2 u = + = + + .
∂x2 ∂y 2 ∂r2 r ∂r r2 ∂θ2
145
El problema que queremos resolver es entonces
∂ 2 u 1 ∂u 1 ∂ 2u
EDP: + + =0 0<r≤R −π ≤θ ≤π
∂r2 r ∂r r2 ∂θ2
CB: u(R, θ) = f (θ), −π ≤ θ ≤ π.
Aún nos falta definir condiciones en θ = ±π. Esto es similar a lo que ocurre en la
situación del alambre circular. El extremo θ = −r corresponde al mismo punto de la
placa circular que θ = π. Aunque allı́ no hay realmente un borde de la placa, el hecho de
que la temperatura sea continua allı́ y que el flujo de calor en la dirección de θ también
sea continuo, implica las siguientes condiciones de borde:
∂u ∂u
u(r, −π) = u(r, π) (r, −π) = (r, π).
∂θ ∂θ
Estas condiciones de borde se conocen como condiciones periódicas. Notamos que estas
condiciones “de borde” son también lineales y homogéneas. De este modo, el problema
luce similar a la ecuación de Laplace en un rectángulo: Hay cuatro condiciones, de las
cuales sólo una es no-homogénea: u(a, θ) = f (θ). Aplicaremos entonces el método de
separación de variables.
Buscaremos soluciones producto
u(r, θ) = φ(θ)G(r),
que satisfagan la EDP y las tres condiciones de borde homogéneas (nuevamente ignoramos
momentáneamente la condición no-homogénea). Las condiciones periódicas implican
Hemos introducido la constante de separación como λ (en lugar de −λ) porque hay dos
condiciones homogéneas para la función φ de la variable θ. El problema que determinará
los autovalores λ es
Los autovalores λ se determinan del modo usual. En efecto, este es uno de los problemas
estándar que ya hemos resuelto, el del alambre circular, con L = π. Por lo tanto los
autovalores son nπ 2
λ= = n2 ,
L
con las autofunciones correspondientes
ó
r2 G00 (r) + rG0 (r) − n2 G(r) = 0. (7.57)
La ecuación (7.57) es lineal y homogénea pero tiene coeficientes no-constantes. Sin
embargo se puede resolver fácilmente. La forma más sencilla de resolver esta ecuación,
es dándose cuenta de que para el operador diferencial en (7.57) la función G(r) = rp se
reproduce a sı́ misma. Esto quiere decir que si reemplazamos rp en el lugar de G(r) en la
ecuación (7.57) obtenemos un múltiplo de rp . En efecto
d2 p d
r2 (r ) + r (rp ) − n2 rp = r2 p (p − 1) rp−2 + r p rp−1 − n2 rp = p (p − 1) + p − n2 rp .
dr 2 dr
Luego nuestro problema se transforma en encontrar las potencias p para las cuales
p(p − 1) + p − n2 = 0.
Pero
p(p − 1) + p − n2 = p2 − p + p − n2 = p2 − n2
y esta última expresión es igual a cero cuando p = n ó p = −n. Si n 6= 0 estas raı́ces son
distintas y las funciones rn y r−n son linealmente independientes, por lo que la solución
general resulta
G(r) = c1 rn + c2 r−n .
147
El caso n = 0 es diferente, e igualmente importante porque λ = 0 también es un autovalor.
Si n = 0 la ecuación (7.57) resulta
r2 G00 (r) + rG0 (r) = 0
que es equivalente a
rG00 (r) + G0 (r) = 0
y también a
d
(rG0 (r)) = 0
dr
lo que a su vez equivale a decir que rG0 (r) es constante. Luego G0 (r) = const/r, integrando
obtenemos
G(r) = c1 + c2 ln r.
Ya hemos discutido la condición en r = 0 y hemos dicho que requeriremos |u(0, θ)| <
∞, lo que implica para la solución producto que
|G(0)| < ∞.
Si ahora observamos las soluciones generales obtenidas, vemos que la parte multiplicada
por c2 tiende a ∞ cuando r → 0. Esto nos dice que c2 = 0 y luego la solución general de
la parte dependiente de r es
G(r) = c1 rn , n ≥ 0,
donde para n = 0 esta expresión se reduce a la constante c1 .
Concluyendo, las soluciones producto que satisfacen las condiciones homogéneas son
rn cos nθ (n ≥ 1) y rn sen nθ (n ≥ 0).
Por el principio de superposición, la solución general de la ecuación de Laplace en el
cı́rculo es ∞ ∞
X X
n
u(r, θ) = A0 + An r cos nθ + Bn rn sen nθ.
n=1 n=1
Este problema es similar al que resolvimos en la Sección 7.3.2. La única diferencia es que
aquı́ el papel de los coeficientes an y bn es jugado por An an y por Bn an , respectivamente.
Usando las mismas fórmulas integrales de esa sección, obtenemos
Z π
1
A0 = f (θ) dθ
2π −π
Z π
1
An an = f (θ) cos nθ dθ, n≥1
π −π
1 π
Z
n
Bn a = f (θ) sen nθ dθ, n ≥ 1.
π −π
148
Como an 6= 0, la solución de la ecuación del calor en el disco de radio a con temperatura
prescripta u(a, θ) = f (θ) en la frontera es
∞
X ∞
X
n
u(r, θ) = A0 + An r cos nθ + Bn rn sen nθ,
n=1 n=1
149
Figura 7.4: Cı́rculo dentro de
una región general
Principio del Máximo. Podemos utilizar este último resultado para probar el princi-
pio del máximo para la ecuación de Laplace: en estado de equilibrio la temperatura
no puede tener un máximo en el interior (a menos que la temperatura sea constante
en toda la región considerada). La demostración se realiza por contradicción. Sea u la
solución de la ecuación de Laplace en una región R. y sea P un punto interior a R, como
el de la Figura 7.4. Por el teorema del valor medio, la temperatura en el punto P es igual
al promedio de la temperatura en el borde del disco de radio r0 . Luego es imposible que
la temperatura en P sea mayor a la temperatura en todos los puntos del borde. Por lo
tanto, u no puede tener un máximo en P . Es decir:
Buen Planteo del Problema de Laplace. El principio del máximo es una herramien-
ta muy importante para el análisis de ecuaciones en derivadas parciales, especialmente al
establecer propiedades cualitativas. Observemos la siguiente definición:
150
Se dice que un problema está bien planteado si existe una única solución que
depende en forma continua de los datos no-homogéneos.
Es decir, un problema estará bien planteado si la solución cambia poco cuando los datos
(valor de borde) cambian poco. Este es un concepto importante en problemas de la Fı́si-
ca. Si la solución cambiara dramáticamente con pequeños cambios en los datos, entonces
cualquier medición fı́sica deberı́a ser exacta para que la solución que se obtiene a partir
de las mediciones sea confiable, o represente bien la realidad. Sabemos que medir exac-
tamente no es posible. Afortunadamente, la mayorı́a de los problemas en ecuaciones en
derivadas parciales están bien planteados.
El principio del máximo puede utilizarse para demostrar que la ecuación de Laplace
2
∇ u = 0 está bien planteada. Supongamos que u y v son dos soluciones de la ecuación
de Laplace con valor de borde f y g, respectivamente, es decir
∇2 u = 0, u = f en el borde,
∇2 v = 0, v = g en el borde.
Supongamos que f y g difieren muy poco y consideremos la diferencia entre estas dos
soluciones: w = u − v. Por la linealidad del operador de Laplace, tenemos que
∇2 w = 0 y w =f −g en el borde.
El principio del máximo (y del mı́nimo) para la ecuación de Laplace implica que el máximo
(y el mı́nimo) de w ocurre en el borde del dominio. Luego, en cualquier punto x interior
se tiene que
mı́n (f (y) − g(y)) ≤ w(x) = u(x) − v(x) ≤ máx (f (y) − g(y)).
y en el borde y en el borde
151
ZZ Z
Por el teorema de la divergencia se tiene que ∇ · (∇u) dx dy = ∇u · n ds, donde
R ∂R
∂R denota la curva cerrada que rodea a R y n el vector normal exterior a ∂R. Luego
Z
∇u · n ds = 0. (7.58)
∂R
Esto nos dice que si u es solución de la ecuación de Laplace (ecuación del calor esta-
cionaria), el flujo de calor neto a través de la frontera debe ser cero. Esto también está
claro desde el punto de vista fı́sico, puesto que en caso contrario habrı́a cambio (en el
tiempo) de la energı́a térmica total dentro de la región, violando la condición de estado
estacionario. La condición (7.58) se denomina condición de compatibilidad para la
ecuación de Laplace.
7.5. Ejercicios
7.1. Consideremos la ecuación
∂u ∂ 2u
= k 2,
∂t ∂x
con las condiciones de borde
u(0, t) = 0 y u(L, t) = 0.
Resolver el problema a valores iniciales si la temperatura es inicialmente
9πx 3πx
(a) u(x, 0) = 6 sen (c) u(x, 0) = 2 cos
L ( L
πx 3πx x si 0 ≤ x ≤ L/2
(b) u(x, 0) = 3 sen − sen (d) u(x, 0) =
L L L − x si L/2 ≤ x ≤ L
Las respuestas a (c) y (d) pueden involucrar ciertas integrales que no es necesario
evaluar, sólo dejarlas expresadas.
En los ejercicios (a) y (b) graficar la temperatura en el punto medio de la barra para
t entre 0 y 10. Determinar el instante de tiempo en que la temperatura en el punto medio
de la barra es la mitad de la temperatura inicial en ese punto.
7.2. Consideremos la ecuación del calor siguiente en una barra de longitud 1:
ut = kuxx ,
0 ≤ x ≤ 1, t ≥ 0,
u(0, t) = 0, u(1, t) = 0, t ≥ 0,
u(x, 0) = 100 sen(πx), 0 ≤ x ≤ 1.
Dar una expresión para la temperatura en el punto medio de la barra en función de
t. ¿Para qué valor de t la temperatura en el punto medio es igual a 50?
7.3. Consideremos la ecuación diferencial
d2 φ
+ λφ = 0.
dx2
Determine los autovalores λ (y autofunciones correspondientes), si φ satisface las siguien-
tes condiciones de borde. Analizar los tres casos (λ > 0, λ = 0, λ < 0).
152
dφ
(a) φ(0) = 0 y φ(π) = 0 (d) φ(0) = 0 y (L) = 0
dx
dφ
(b) φ(0) = 0 y φ(1) = 0 (e) (0) = 0 y φ(L) = 0
dx
dφ dφ dφ dφ
(c) (0) = 0 y (L) = 0 (f ) φ(−L) = φ(L) y (−L) = (L)
dx dx dx dx
∂u ∂ 2u
7.4. Resolver la ecuación del calor = k 2 , 0 ≤ x ≤ L, t ≥ 0, con las condiciones de
∂t ∂x
borde
∂u ∂u
(0, t) = 0, (L, t) = 0, t > 0,
∂x ∂x
y con la condición inicial siguiente:
3πx πx
(a) u(x, 0) = 6 + 4 cos (c) u(x, 0) = −2 sen
L ( L
8πx 0 si x < L/2
(b) u(x, 0) = −3 cos (d) u(x, 0) =
L 1 si x > L/2
En los ejercicios (a) y (b) graficar la temperatura en el extremo derecho de la barra
para t entre 0 y 10. Determinar el instante de tiempo en que la temperatura en el extremo
derecho de la barra es la mitad de la temperatura inicial en ese punto. En cada inciso,
determinar la temperatura cuando t → ∞.
7.5. Consideremos la ecuación del calor siguiente en una barra de longitud 1:
1
ut = 10 uxx ,
0 ≤ x ≤ 1, t ≥ 0,
ux (0, t) = 0, ux (1, t) = 0, t ≥ 0,
u(x, 0) = f (x), 0 ≤ x ≤ 1.
7.6. Resolver
∂u ∂ 2u ∂u
=k 2 con (0, t) = 0
∂t ∂x ∂x
u(L, t) = 0
u(x, 0) = f (x).
153
7.7. Resolver la ecuación de Laplace dentro de un rectángulo 0 ≤ x ≤ L, 0 ≤ y ≤ H con
las siguientes condiciones de borde:
∂u ∂u
(a) (0, y) = 0, (L, y) = 0, u(x, 0) = 0, u(x, H) = f (x).
∂x ∂x
∂u ∂u
(b) (0, y) = g(y), (L, y) = 0, u(x, 0) = 0, u(x, H) = 0.
∂x ∂x
∂u
(c) (0, y) = 0, u(L, y) = g(y), u(x, 0) = 0, u(x, H) = 0.
∂x
∂u
(d) u(0, y) = g(y), u(L, y) = 0, (x, 0) = 0, u(x, H) = 0.
∂y
∂u ∂u
(e) u(0, y) = f (y), u(L, y) = 0, (x, 0) = 0, (x, H) = 0.
∂y ∂y
(
∂u ∂u 0 si x > L/2, ∂u
(f ) (0, y) = 0, (L, y) = 0, u(x, 0) = (x, H) = 0.
∂x ∂x 1 si x < L/2, ∂y
7.8. Supongamos que u(x, y) es la solución de la ecuación de Laplace dentro del rectángulo
0 ≤ x ≤ L, 0 ≤ y ≤ H con las siguientes condiciones de borde:
∂u ∂u ∂u ∂u
(0, y) = 0 (x, 0) = 0 (L, y) = 0 (x, H) = f (x).
∂x ∂y ∂x ∂y
(a) Sin resolver el problema, explicar la condición que debe cumplir f para que este
problema tenga solución.
(b) Resolver el problema por el método de separación de variables. Notar que queda una
constante sin determinar, es decir que puede tomar cualquier valor.
7.9. Resolver la ecuación de Laplace fuera de un cı́rculo de radio a con la siguiente
condición de borde (suponer que u(r, θ) permanece acotado cuando r → ∞.)
(a) u(a, θ) = ln 2 + 4 cos 3θ
154
(a) u = 0 sobre el diámetro (parte recta) y u(a, θ) = g(θ)
(b) el diámetro está aislado y u(a, θ) = g(θ).
7.12. Resolver la ecuación de Laplace en un anillo circular a ≤ r ≤ b con las siguientes
condiciones de borde
(a) u(a, θ) = f (θ), u(b, θ) = 0
∂u
(b) (a, θ) = 0, u(b, θ) = f (θ)
∂r
7.13. Resolver la ecuación de Laplace en un sector de 90◦ de un anillo circular (a ≤ r ≤ b,
0 ≤ θ ≤ π/2) con las siguientes condiciones de borde:
u(r, 0) = 0, u(r, π/2) = 0, u(a, θ) = 0, u(b, θ) = f (θ).
p
7.14. Consideremos el siguiente problema, planteado sobre la región anular 1 ≤ x2 + y 2 ≤
2:
∇2 u = 0 1 ≤ r ≤ 2,
u=0 sobre el borde r = 1,
u = 100 sobre el borde r = 2.
(a) Hallar la solución (ayuda, es de variables separables en coordenadas polares r, θ).
(b) Determinar la curva sobre la que u = 50.
7.15. Consideremos la ecuación de Laplace siguiente en el cuadrado [0, 1] × [0, 1]:
uxx + uyy = 0, 0 ≤ x ≤ 1, 0 ≤ y ≤ 1,
0 ≤ y ≤ 1,
ux (0, y) = 0, u(1, y) = 0,
u(x, 0) = 0, u(x, 1) = 30 cos( πx ) − 10 cos( 3πx ),
0 ≤ x ≤ 1.
2 2
(a) Halle la solución u(x, y).
(b) Dé una fórmula para la solución u sobre la lı́nea x = 0.2, (0 ≤ y ≤ 1). Grafı́quela en
la computadora.
(c) Determine el valor aproximado (5 dı́gitos) de y para el que u(0.2, y) = 15.
7.16. Consideremos la ecuación de Laplace siguiente en el cuadrado [0, 1] × [0, 1]:
uxx + uyy = 0, 0 ≤ x ≤ 1, 0 ≤ y ≤ 1,
0 ≤ y ≤ 1,
u(0, y) = 0, ux (1, y) = 0,
πx 3πx
u(x, 0) = 0, u(x, 1) = 30 sen( ) − 10 sen( ), 0 ≤ x ≤ 1.
2 2
(a) Diga cuál es la solución u(x, y).
(b) Dé una fórmula para la solución u sobre la lı́nea x = 0.9, (0 ≤ y ≤ 1). Grafı́quela en
la computadora.
(c) Determine el valor aproximado (5 dı́gitos) de y para el que u(0.9, y) = 30.
155
Capı́tulo 8
Diferencias Finitas
156
Idea: La idea principal del método de diferencias finitas consiste en reemplazar
derivadas por cocientes de diferencias.
dg g(x + h) − g(x)
Recordando que (x) = lı́m , la idea es tomar h pequeño y aproxi-
dx h→0 h
dg
mar dx (x) por g(x+h)−g(x)
h
, o alguna variante:
dg
Diferencia Fórmula Error: fórmula −
dx
dg g(x + h) − g(x)
Adelantada: (x) ≈ O(h) si g ∈ C 2
dx h
dg g(x) − g(x − h)
Atrasada: (x) ≈ O(h) si g ∈ C 2
dx h
dg g(x + h2 ) − g(x − h2 )
Centrada: (x) ≈ O(h2 ) si g ∈ C 3
dx h
Las acotaciones del error se obtienen utilizando el Teorema de Taylor. Recordemos
2
que en la ecuación de Poisson aparece uxx = ddxu2 . Para derivadas segundas hacemos lo
siguiente, utilizando diferencias centradas:
du
d2 u + h2 ) − du
− h2 )
d du dx
(x dx
(x
2
(x) = (x) ≈
dx dx dx h
u(x+ h +h )−u(x+ h −h ) u(x− h +h )−u(x− h −h )
2 2
h
2 2
− 2 2
h
2 2
≈
h
u(x+h)−u(x)
h
− u(x)−u(x−h)
h
=
h
u(x + h) − 2u(x) + u(x − h)
=
h2
Para esta aproximación de la derivada segunda tenemos el siguiente resultado:
h2 00 h3 h4
u(x + h) = u(x) + hu0 (x) + u (x) + u000 (x) + uIV (ξ1 )
2 3! 4!
2 3 4
h h h
u(x − h) = u(x) − hu0 (x) + u00 (x) − u000 (x) + uIV (ξ2 ),
2 3! 4!
157
con x − h < ξ2 < x < ξ1 < x + h. Luego, sumando ambas igualdades, obtenemos
u(x1 ) = a,
u(xi+1 ) − 2u(xi ) + u(xi−1 )
−k ≈ f (xi ), i = 2, 3, . . . , N,
h2
UN +1 = b.
158
Este sistema de ecuaciones tiene N + 1 ecuaciones y N + 1 incógnitas U1 , U2 , . . . , UN +1 .
El mismo, puede reescribirse de la siguiente manera:
U1 =a
h2
−U + 2U − U = f2
1 2 3
k
h2
− U2 + 2U3 − U4 = f3
k
. .. .
..
h2
−U + 2U − U = fN
N −1 N N +1
k
UN +1 = b,
Script poissondirichletdf1d.m.
% poissondirichletdf1d.m
%
% Programita para resolver la ecuacion de Poisson
% con condiciones de Dirichlet.
% Metodo: Diferencias finitas
%
% - k u’’(x) = f(x), 0 < x < L
% u(0) = a, u(L) = b
% % Armado de la matriz
h = L/N;
unos = ones(N+1,1);
diagonales = [-1*unos 2*unos -1*unos];
159
matriz = spdiags(diagonales, [-1 0 1], N+1, N+1);
matriz(1,[1:2]) = [1 0 ];
matriz(N+1,[N:N+1]) = [0 1];
% % Armado del lado derecho
X = linspace(0,L,N+1)’;
F = [a ; h^2/k*f(X(2:N)) ; b];
% % Resolucion
% Primero resolvemos para los puntos interiores (i=2:N)
U = matriz \ F;
% Ahora agregamos los valores en los extremos para graficar
figure(1); plot(X,U,’*-’)
Lo que puede demostrarse, a partir del error cometido por la fórmula en diferen-
cias (8.2) es que
Teorema 56. Si la solución exacta u(x) de (8.1) es C 4 [0, L] entonces
1
máx |u(xi ) − Ui | ≤ CM4 h2 = CM4 L2 ,
i=1,2,...,N +1 N2
donde M4 = máx[0,L] |uIV | y C es una constante que depende de los parámetros k, L, de
la ecuación, pero es independiente de la función f (x), y de los datos de borde a, b y del
parámetro de discretización h (o N ).
Demostración. La demostración de este teorema queda fuera del alcance de este curso.
El lector interesado puede encontrarla en [Larsson-Thomée 2009, Ch. 4]
Observación 57. Cuando uno escribe un programa como el anterior para resolver un
problema, debe verificar que funcione bien. Para ello, se busca una solución exacta. Esto
parece difı́cil, pero es cuestión de elegir una u(x) cualquiera, que sea C 4 [0, L] aunque no
sea muy obvia, como lineal o cuadrática y luego, dado un k elegido, se calcula f (x), u(0)
y u(L).
2
Por ejemplo: Tomemos L = 1, k = 2 y u(x) = x2 + e−(x−0.5) . Para que esta u sea
solución de (8.1), calculamos f = −kuxx :
2 2
ux (x) = 2x − 2(x − 0.5)e−(x−0.5) , uxx (x) = 2 − 2[1 − 2(x − 0.5)2 ]e−(x−0.5)
Luego, resolvemos el problema con
2 −(x−0.5)2
f (x) = −4 + 4[1 − 2(x − 0.5) ]e = −kuxx (x) con k = 2
2
a = e−0.5 = u(0)
2
b = 1 + e−0.5 = u(1) con L = 1
2
y deberı́amos obtener una solución aproximada a u(x) = x2 + e−(x−0.5) .
Para estar seguros que el código está bien programado, debemos verificar que cuando
tomamos h más pequeño, el error máximo se reduce como Ch2 . Es decir, cada vez que
tomamos h igual a la mitad de un h anterior, el error debe reducirse por un factor 1/4.
160
Otras Condiciones de Borde
En esta sección consideramos en el extremo derecho x = L una condición de borde
diferente, que contiene a la de tipo Neumann y de tipo Robin.
Consideremos el problema
−kuxx = f, 0 ≤ x ≤ L,
u(0) = a, (8.3)
ku0 (L) + H u(L) = H u .
1 2 E
u(L + h) − u(L − h)
u0 (L) ≈ .
2h
Notamos que xN +2 = L + h es un punto que cae fuera del dominio [0, L] y por eso a xN +2
lo llamamos nodo ficticio. La condición de borde se reemplaza entonces por la ecuación
UN +2 − UN
k + H1 UN +1 = H2 uE .
2h
Esta ecuación reemplaza a la ecuación UN +1 = b en el sistema del problema anterior.
Observamos que seguimos teniendo N +1 ecuaciones, pero ahora tenemos N +2 incógnitas,
por lo que nos falta una ecuación más. Lo que hacemos es considerar xN +1 = L como un
nodo interior : imponemos la ecuación diferencial −ku00 = f en x = xN +1 = L. Es decir,
agregamos al sistema la ecuación
−UN + 2UN +1 − UN +2
k = fN +1 .
h2
De esta manera, el sistema resulta
U1 =a
h2
−U1 + 2U2 − U3
= f2
k
h2
− U2 + 2U3 − U4 = f3
k
.. ..
. . (8.4)
2
h
−UN −1 + 2UN − UN +1
= fN
k
h2
−UN + 2UN +1 − UN +2 = fN +1
k
2hH1 2hH2
−UN + UN +1 + UN +2 = uE ,
k k
161
De forma matricial se escribe de la siguiente manera:
a
1 0 0 ... ... 0
U1
h2
−1 U2 k 2
f
−1 0
2 ... 0
h2
.. .. U3
k 3
f
0
−1 2 −1 . .
..
..
= .
.
.. .. .. .. .. . .
. . . . 0 U h2
N 2k f N
0 . . . 0 −1 2 −1 UN +1 fN +1
h
2hH1
0 . . . 0 −1 k
1 UN +2
k
2hH2
uE k
Script poissonrobindf1d.m.
% poissonrobindf1d.m
%
% Programita para resolver la ecuacion de Poisson
% con condiciones de Dirichlet en x=0 y Robin en x=L.
% Metodo: Diferencias finitas
%
% - k u’’(x) = f(x), 0 < x < L
% u(0) = a, k u’(L) + H1 u(L) = H2 uE
% % Armado de la matriz
h = L/N;
unos = ones(N+2,1);
diagonales = [-1*unos 2*unos -1*unos];
matriz = spdiags(diagonales, [-1 0 1], N+2, N+2);
matriz(1,[1:2]) = [1 0];
matriz(N+2, N:N+2) = [-1 2*h*H1/k 1];
% % Resolucion
U = matriz \ F;
162
% Ahora eliminamos el valor del punto ficticio N+2
U(N+2) = []; % otra forma: U = U(1:N+1);
figure(1); plot(X,U,’m-’)
También para este problema se tiene un resultado de estimación del error, cuya de-
mostración escapa al alcance de este curso:
Aquı́ los datos del problema son: k, L, T , f (x, t), a(t), H1 , H2 , uE (t) y u0 (x).
Para resolver este problema consideramos la grilla de puntos (xi , tj ) con xi = (i − 1)h
como antes y tj = j∆t, donde ∆t es el parámetro de discretización temporal.
Método explı́cito
Si aproximamos ut con una diferencia adelantada
Esta fórmula será útil para calcular U en el interior de [0, L], para t > 0.
En t = 0 (j = 0) usaremos la condición inicial u(x, 0) = u0 (x), por lo que tomaremos
0
Ui = u0 (xi ), i = 1, 2, . . . , N + 1.
En x = 0 usaremos la condición de borde u(0, t) = a(t) por lo que tomaremos U1j =
a(tj ).
En x = L tenemos la condición de borde kux (L, t) + H1 u(L, t) = H2 uE (t). Aproxi-
mando ux (L, t) ≈ u(L+h,t)−u(L−h,t)
h
llegamos a la fórmula
UNj +2 − UNj 2h
k + H1 UNj +1 = H2 uE (tj ) UNj +2 = (H2 uE (tj ) − H1 UNj +1 ) + UNj .
2h k
A partir de estas fórmula, diseñamos el siguiente método numérico para la ecuación
de difusión (8.5): Sean N, M ∈ N y definamos h = L/N , ∆t = T /M , λ = k∆t/h2
(j = 0) Ui0 = u0 (xi ), i = 1, 2, . . . , N + 1;
Para j = 0, 2, . . . M − 1
• U1j+1 = a(tj+1 );
• Uij+1 = ∆tfij + λUi−1
j
+ (1 − 2λ)Uij + λUi+1
j
, i = 2, . . . , N + 1;
• UNj+1
+2 =
2h
k
(H2 uE (tj+1 ) − H1 UNj+1 j+1
+1 ) + UN ;
Script difusionrobindf1dexpl.m.
% difusion_robin_df_1d_expl.m
%
% Programita para resolver la ecuacion de Difusion
% con condiciones de Dirichlet en x=0 y Robin en x=L.
% Metodo: Diferencias finitas en espacio y Euler EXPLICITO en tiempo.
%
% u_t(x,t) - k u_xx(x,t) = f(x,t), 0 < x < L, t > 0
% u(0,t) = a(t), k u’(L,t) + H1 u(L,t) = H2 uE(t), t > 0
% u(x,0) = u0(x), 0 < x < L.
164
% % Parametros del problema
L = 1; T = 2.5; k = 1;
f = @(x,t)(10*exp(-100*(x-.5).^2)*(t<1.5));
a = @(t)(0);
H1 = 3; H2 = 3; uE = @(t)(0);
u0 = @(x)(zeros(size(x)));
% % Condicion inicial
X = linspace(0,L,N+1)’;
U = u0(X);
% le damos valor a la U en el nodo ficticio (t=0)
U(N+2) = (H2*uE(0)-H1*U(N+1))*2*h/k + U(N);
ejes = [0 L 0 1];
figure(1); plot(X,U(1:N+1));
title(’t = 0’)
axis(ejes); pause
% % Resolucion
for t=deltat:deltat:T
% calculamos U en los nodos "interiores"
U(2:N+1) = deltat*f(X(2:N+1), t) ...
+ lambda*U(1:N) + (1-2*lambda)*U(2:N+1) + lambda*U(3:N+2);
% le damos valor a U en el extremo izquierdo
U(1) = a(t);
% le damos valor a la U en el nodo ficticio
U(N+2) = (H2*uE(t)-H1*U(N+1))*2*h/k + U(N);
figure(1); plot(X,U(1:N+1),’*-’);
title(sprintf(’t = %5.3f’,t));
axis(ejes); pause(0.01);
end
Este método tiene una restricción para resultar estable. Si observamos el término que
multiplica Uij en la fórmula para Uij+1 vemos que es 1 − 2λ con λ = k∆t/h2 . Resulta que
si λ ≥ 1/2, de manera que 1 − 2λ ≤ 0 se observa un fenómeno de inestabilidad (verificarlo
con el código) que hace que las soluciones oscilen en cada paso de tiempo. Lo que ocurre
es que pequeños errores de redondeo se magnifican exponencialmente.
Se recomienda, una vez elegido h, elegir ∆t de manera que λ ≤ 1/4, para que el
algoritmo funcione correctamente. Más aún, en el código que mostramos, elegimos λ y
165
2
luego ∆t se calcula para que λ = k∆t
h2
, es decir ∆t = λ hk .
También hay estimaciones del error, que resumimos en el siguiente Teorema.
Teorema 59. Si la solución exacta u(x, t) de (8.5) es C 4 [0, L] para la variable x y C 2 [0, T ]
para la variable T , entonces, si λ = k∆t/h2 < 1/2, resulta
donde M4,2 = máx[0,L]×[0,T ] |uxxxx |+|utt | y C es una constante que depende de los paráme-
tros k, L, H1 , H2 de la ecuación, pero es independiente de la función f (x, t), del dato
inicial u0 (x), de los datos de borde a(t), uE (t) y de los parámetro de discretización h, ∆t.
Método implı́cito
Cuando los coeficientes son variables, o cuando queremos resolver problemas en dos
o tres dimensiones espaciales, la condición de estabilidad puede ser más difı́cil de deter-
minar. En esta sección veremos un método que es incondicionalmente estable, es decir,
resulta estable sin importar la relación entre ∆t y h.
La diferencia principal con el método obtenido anteriormente está en aproximar la
derivada temporal ut por una diferencia atrasada en lugar de adelantada. Es decir, con-
sideraremos que
u(xi , tj+1 ) − u(xi , tj )
≈ ut (xi , tj+1 ).
∆t
Por lo tanto, vemos que
Vemos que no queda Uij+1 despejada sola a la izquierda, sino que queda relacionada con
j+1 j+1
Ui−1 y con Ui+1 . Por eso se dice que el método es implı́cito.
Por lo tanto, en cada paso de tiempo deberemos resolver un sistema de ecuaciones.
Esto es un poco más costoso computacionalmente que lo que se hace en el caso explı́cito,
pero se gana en estabilidad.
Si incorporamos las condiciones de borde, vemos que, conocido Uij , para un j dado y
para i = 1, 2, . . . , N + 2, los valores Uij+1 , i = 1, 2, . . . , N + 2 deben satisfacer el siguiente
166
sistema de N + 2 ecuaciones
U1j+1 = a(tj+1 ),
j+1
+ (1 + 2λ)Uij+1 − λUi+1
j+1
= ∆tfij+1 + Uij ,
−λUi−1 i = 2, 3, . . . , N + 1, (8.6)
2hH1 j+1 2hH2
−UNj+1 + UN +1 + UNj+1
+2 = uE (tj+1 ).
k k
En forma matricial, el sistema a resolver en cada paso de tiempo resulta:
U1j+1
a(tj+1 ) 0
1 0 0 ... ... 0
j+1 j f2j+1
−λ 1 + 2λ −λ U2j+1 U2j
0 . . . 0 j+1
U3 U3
.. f
...
3
0 −λ 1 + 2λ −λ . .. ..
.
. . . . . . = . + ∆t .. .
.. .. .. .. .. 0 UNj+1 UNj
j+1
fN
0 ... 0 −λ 1 + 2λ −λ j+1 j
j+1
U N +1 U N +1
f
N +1
0 ... 0 −1 2hH 1
1
k UNj+1
+2
2hH2
k
uE 0
Script difusionrobindf1dimpl.m.
% difusion_robin_df_impl_1d.m
%
% Programita para resolver la ecuacion de Difusion
% con condiciones de Dirichlet en x=0 y Robin en x=L.
% Metodo: Diferencias finitas en espacio y Euler IMPLICITO en tiempo.
%
% u_t(x,t) - k u_xx(x,t) = f(x,t), 0 < x < L, t > 0,
% u(0,t) = a(t), k u’(L,t) + H1 u(L,t) = H2 uE(t), t > 0,
% u(x,0) = u0(x), 0 < x < L.
167
% % Condicion inicial
X = linspace(0,L,N+1)’;
U = u0(X);
% le damos valor a la U en el nodo ficticio
U(N+2) = (H2*uE(0)-H1*U(N+1))*2*h/k + U(N);
ejes = [0 L 0 1];
figure(1); plot(X,U(1:N+1));
title(’t = 0’);
axis(ejes); pause
% % Armado de la matriz
unos = ones(N+2,1);
columnas = [-lambda*unos (1+2*lambda)*unos -lambda*unos];
% % Resolucion
F=zeros(N+2,1);
for t=deltat:deltat:T
% Ensamblado del lado derecho
U(1) = a(t);
U(N+2) = 2*H2*h/k*uE(t);
F = [ 0 ; f(X(2:N+1),t) ; 0 ];
lado_derecho = U + deltat*F;
% calculamos la nueva U
U = matriz \ lado_derecho;
figure(1);
plot(X,U(1:N+1),’*-’);
title(sprintf(’t = %5.3f’,t));
axis(ejes); pause(0.01);
end
Teorema 60. Si la solución exacta u(x, t) de (8.5) es C 4 [0, L] para la variable x y C 2 [0, T ]
168
para la variable T , entonces resulta
donde M = máx[0,L]×[0,T ] |uxxxx |+|utt | y C es una constante que depende de los parámetros
k, L, H1 , H2 de la ecuación, pero es independiente de la función f (x, t), del dato inicial
u0 (x), de los datos de borde a(t), uE (t) y de los parámetro de discretización h, ∆t.
∂ 2u ∂ 2u
∇2 u(xi , yj ) = (x i , yj ) + (xi , yj )
∂x2 ∂y 2
ui+1,j − 2ui,j + ui−1,j h2 ∂ 4 u ui,j+1 − 2ui,j + ui,j−1 h2 ∂ 4 u
= + (ξi,j , yj ) + + (xi , ηi,j )
h2 12 ∂x4 h2 12 ∂y 4
ui+1,j + ui,j+1 − 4ui,j + ui−1,j + ui,j−1 h2 h ∂ 4 u ∂ 4u i
= + (ξi,j , yj ) + (x i , ηi,j )
h2 12 ∂x4 ∂y 4
Por lo tanto, si la solución u(x, y) de (8.7) es C 4 en Ω̄, satisface
−ui+1,j − ui,j+1 + 4ui,j − ui−1,j − ui,j−1
k = fi,j + O(h2 )
h2
para i, j = 1, 2, . . . , N −1. Aquı́ fi,j denota la evaluación de f en (xi , yj ), respectivamente.
169
A partir de estas expresiones, definimos Ui,j , aproximaciones de ui,j a través de las
siguientes igualdades:
−Ui+1,j − Ui,j+1 + 4Ui,j − Ui−1,j − Ui,j−1
k = gi,j , i, j = 1, 2, . . . , N − 1,
h2
(8.8)
U0,j = α(0, yj ), UN,j = α(1, yj ), j = 0, 1, . . . , N,
Ui,0 = α(xi , 0), Ui,N = α(xi , 1), i = 0, 1, . . . , N.
Observación 61. Para implementarlo debemos numerar los nodos e incógnitas con un
solo ı́ndice, por ejemplo k = (i − 1)(N + 1) + j. Al hacer esto, el sistema resulta con
una matriz pentadiagonal en lugar de tridiagonal. A continuación mostramos un código
ejemplo
% poisson_dirichlet_2d.m
% Metodo de diferencias finitas en 2d
% para el problema
%
% - k (u_xx + u_yy) = f(x,y), (x,y) en Omega = (a,b) x (a,b)
% u(x,y) = alfa(x,y), (x,y) en la frontera de Omega
%
%% Parametros de la resolucion
N = 30;
%% Resolucion
h=1/N;
xx=linspace(a,b,N+1);
yy=xx;
[x,y] = meshgrid(xx,yy);
x = x’; x=x(:);
y = y’; y=y(:);
F=f(x,y);
ALFA=alfa(x,y);
170
M = spdiags(diagonales, [N+1 1 0 -1 -(N+1)], (N+1)^2, (N+1)^2);
M = M’;
rhs = F;
%% Resolvemos
U = M\rhs;
171
otro problema del mismo tipo, y ası́ sucesivamente. Cada uno de estos problemas puede
resolverse de manera similar a la ecuación de Poisson, discretizando en x y en y con un
parámetro h = 1/N . Más precisamente, aproximaremos unxx (xi , yj ) + unyy (xi , yj ) por la
fórmula
uni+1,j + uni,j+1 − 4uni,j + uni−1,j + uni,j−1
h2
Reemplazando esta fórmula en (8.9) obtenemos el siguiente sistema de ecuaciones para
n
Ui,j , que son las aproximaciones de uni,j , para n = 1, 2, 3, . . .
n k∆t n n n n n
n−1 n
Ui,j − 2
U i+1,j + U i,j+1 − 4Ui,j + U i−1,j + Ui,j−1 = Ui,j + ∆tfi,j , i, j = 1, 2, . . . , N − 1,
h
n
U0,j = α(0, yj , tn ), UN,j = α(1, yj , tn ), j = 0, 1, . . . , N,
Ui,0 = α(xi , 0, tn ), Ui,N = α(xi , 1, tn ), i = 0, 1, . . . , N.
k∆t
Llamando λ = h2
obtenemos
n n n n n n−1 n
− λUi+1,j − λUi,j+1 + (1 + 4λ)4Ui,j − λUi−1,j − λUi,j−1 = Ui,j + ∆tfi,j , i, j = 1, 2, . . . , N − 1,
n n n
U0,j = α(0, yj , t ), UN,j = α(1, yj , t ), j = 0, 1, . . . , N,
Ui,0 = α(xi , 0, tn ), Ui,N = α(xi , 1, tn ), i = 0, 1, . . . , N.
Esto resulta en un sistema de (N + 1)2 ecuaciones con (N + 1)2 incógnitas, que debe ser
resuelto en cada paso de tiempo. Es decir, para n = 1, 2, . . . .
Sin decirlo explı́citamente, hemos supuesto que la condición inicial viene dada por
0
Ui,j = u0 (xi , yj ). La implementación de este método se puede encontrar en el entorno
virtual, en el archivo: difusion_dirichlet_df_2d.m
8.3. Ejercicios
8.1. Demuestre que si u ∈ C 2 , entonces
0
u (x) − u(x + h) − u(x) M2
≤ h,
h 2
172
8.4. Consideremos el siguiente problema:
−10(x−0.7)2
−uxx = 20e
, 0 ≤ x ≤ 1,
u(0) = 5,
u(1) = 6.
Grafique las soluciones obtenidas con N = 100 (h = 1/100). ¿Se anima a dar una interpre-
tación fı́sica de las soluciones obtenidas? Recuerde que la ecuación estudiada corresponde
al estado estacionario de la ecuación de difusión.
8.6. Escriba una discretización análoga a (8.4) para el problema estacionario de difusión
y reacción siguiente:
−kuxx + c u = f, 0 ≤ x ≤ L,
u(0) = a,
ku0 (L) + H u(L) = H u .
1 2 E
Calcule utilizando la computadora y dé una aproximación de u(1/2) con un error absoluto
menor a 10−3 .
173
8.9. Consideremos la siguiente ecuación del calor:
ut − uxx = 100(1 − x), 0 ≤ x ≤ 1, t ≥ 0,
u(0, t) = 0, t ≥ 0,
ux (1, t) = 0, t ≥ 0,
u(x, 0) = 5x(2 − x), 0 ≤ x ≤ 1.
Determine en cada caso la solución estacionaria u∞ (x) y el tiempo t para el cual se cumple
A partir de estos cálculos, deduzca qué efecto tiene sobre la solución el cambio en el
coeficiente de difusión.
8.11. Modifique el programa difusionrobindf1dimpl.m para resolver problemas de
reacción-difusión como el siguiente:
ut − kuxx + c u = f, 0 ≤ x ≤ L, t ≥ 0,
u(0, t) = a(t) t ≥ 0,
kux (L, t) + H1 u(L, t) = H2 uE (t), t ≥ 0,
u(x, 0) = u0 (x), 0 ≤ x ≤ L,
174
8.13. Considere la siguiente ecuación del calor:
ut − uxx = 100, 0 ≤ x ≤ 1, t ≥ 0,
u(0, t) = 0, t ≥ 0,
(∗)
ux (1, t) = −u(1, t), t ≥ 0,
u(x, 0) = 10x(1 − x)2 , 0 ≤ x ≤ 1.
(b) Para u(x, t) la solución de (∗), determine el lı́mite cuando t → ∞ de u(1/2, t).
(a) Modifique
( el código desarrollado en clase para resolver este problema con f (x, t) =
100, si t < 2,
. Resuelva con N = 40 y λ = 4. Visualice la solución durante el
0, si t ≥ 2.
intervalo de tiempo entre t = 0 y t = 4. Explique con palabras lo que se observa,
interpretando u como una temperatura.
Si le interesa leer un poco más sobre los temas de este capı́tulo, le recomendamos el
siguiente libro:
175
Bibliografı́a
[Larsson-Thomée 2009] Larsson, S., Thomée, V., Partial Differential Equations with Nu-
merical Methods, Springer, 2009.
176