PDF de programación - Integración numérica y resolución de algunas ecuaciones diferenciales ordinarias con Maxima CAS

Integración numérica y resolución de algunas ecuaciones diferenciales

ordinarias con Maxima CAS

Renato Álvarez Nodarse

Dpto. Análisis Matemático, Universidad de Sevilla.

26 de junio de 2017

http://euler.us.es/~renato/

Fórmulas de cuadratura.

b

a

Sea f (x) una función continua definida en el intervalo [a, b]. Nuestro objetivo será encontrar

b

fórmulas aproximadas para calcular la integral

f (x) dx. En caso de conocer la primitiva F (x) es

evidente que podemos encontrar el valor exacto de la integral utilizando el Teorema fundamental del
f (x) dx = F (b) − F (a). Sin embargo no siempre esto es posible. Por ejemplo,
cálculo integral:
para la función f (x) = e−x2 no existe ninguna primitiva que podamos escribir utilizando funciones
elementales. En esta práctica vamos a aprender tres métodos para calcular aproximadamente el valor
númerico de las integrales definidas.

a

Fórmula de los rectángulos.

Una aproximación de la integral

por un rectángulo de base b − a y altura f a+b

a

b

b

a

f (x) dx = (b − a)f

(ver figura 1), entonces
a + b



2

+ R(ξ),

2

ξ ∈ [a, b],

f (x) dx consiste en aproximar el área bajo la curva y = f (x)

(1)

donde el error R(ξ), si f tiene primera y segunda derivadas continuas en [a, b], se expresa de la forma

R(ξ) =

f(ξ),

ξ ∈ [a, b].

(2)

(b − a)2

24

b

Ahora si queremos aproximar la integral

f (x) dx con mejor exactitud podemos dividir el inter-

valo [a, b] en n puntos, o sea, consideremos la partición del intervalo

a

[a, b] = [a, x1] ∪ [x1, x2] ∪ ··· ∪ [xn−2, xn−1] ∪ [xn−1, b],

donde

De

b

a

xk = a +

k, n = 0, 1, 2, ..., n,

x0 = a, xn = b.

f (x) dx =

f (x) dx + ··· +

f (x) dx + ··· +

b − a
n

x1

xk+1

a

xk

1

b

f (x) dx.

xn−1

Figura 1: Aproximación de una integral por el método de los rectángulos. A la izquierda vemos el
área bajo la curva que queremos calcular. A la derecha, la aproximación mediante el correspondiente
rectángulo.

si aplicamos a cada integral

f (x) dx la fórmula (1)obtenemos la ecuación

xk+1
b

xk

n−1

k=0

f

b − a
n

xk + xk+1



2

f (x) dx =

+ R(ξ),

|R(ξ)| ≤ M

(b − a)2
24n2

, M = máx
x∈[a,b]

|f(x)|.

y

a

(3)

(4)

Problema 1 Utilizando las fórmulas (1) y (2) demostrar las fórmulas (3) y (4).

Problema 2
Prueba la fórmula (1) y (2). Para ello usa el teorema de Taylor para aproximar f (x)
2 así como el teorema del valor medio integral: Si f : [a, b] −→ R es continua en [a, b]
en el punto a+b
y g : [a, b] −→ R es una función no positiva (no negativa) e integrable en [a, b] entonces existe un
ξ ∈ [a, b] tal que

b

b

f (x)g(x) dx = f (ξ)

g(x) dx.

(5)

Implementación con Maxima.

a

a

Veamos como podemos implementar lo anterior con Maxima CAS. Definimos el intervalo [a, b],

el número de puntos en que vamos a dividir en intervalo y definimos la partición que usaremos:

a:0;b:1;

(%i1)
x[0]:a;
Nu:20;
x[n]:=x[0]+n*(b-a)/Nu;
(%o1)
(%o2)
(%o3)
(%o4)
(%o5)

0
1
0
20
x[n]:=x[0]+(n*(b−a))/Nu

(%i6) define(f(x),x^2);
rec:sum(f((x[k]+x[k+1])/2),k,0,Nu-1)*((b-a)/Nu);
float(%);
(%o6)
(%o7)
(%o8)

f(x):=x^2
533/1600
0.333125

2

yxf(x)yxf(x)En este caso, como la función tiene primitiva podemos comparar el resultado numérico con el valor
exacto

(%i9) exac:float(integrate(f(x),x,a,b));
float(abs(rec-exac));
(%o9)
(%o10) 2.083333333333104*10^−4

0.33333333333333

Fórmula de los trapecios.

Otra aproximación de la integral

f (x) dx consiste en aproximar el área bajo la curva y = f (x)

no por un rectángulo sino por un trapecio de base b − a (ver figura 2), entonces

a

f (x) dx = (b − a)

+ R(ξ),

(6)

b

b

a

f (a) + f (b)



2

donde el error R(ξ), si f tiene primera y segunda derivadas continuas en [a, b] se expresa de la forma

R(ξ) = − (b − a)2

12

f(ξ),

ξ ∈ [a, b].

(7)

Figura 2: Aproximación de una integral por el método de los trapecios. A la izquierda vemos el
área bajo la curva que queremos calcular. A la derecha, la aproximación mediante el correspondiente
trapecio.

Problema 3 Demostrar las fórmulas (6) y (7). Para ello seguir los siguientes pasos:

1. Demostrar que b

a

f(x)(x − a)(x − b) dx = −(b − a)[f (a) + f (b)] + 2

b

a

f (x) dx.

(8)

2. Utilizando el teorema del valor medio integral (5) demostrar que

f(x)(x − a)(x − b) dx = − (b − a)3

6

f(ξ),

ξ ∈ [a, b].

(9)

b

a

3. Usando los dos apartados anteriores obtén las fórmulas (6) y (7).

b

Ahora podemos aproximar la integral

f (x) dx con mejor exactitud dividiendo, igual que antes,

el intervalo [a, b] en n puntos, o sea, consideremos la partición del intervalo
[a, b] = [a, x1] ∪ [x1, x2] ∪ ··· ∪ [xn−2, xn−1] ∪ [xn−1, b],

a

3

yxf(x)yxf(x)donde

Nuevamente, b

b − a
n

k,

x1

b

f (x) dx,

xn−1

y, por tanto, si aplicamos a cada integral

f (x) dx la fórmula (1) obtenemos la expresión

xk = a +

k = 0, 1, 2, ..., n,

x0 = a, xn = b.

xk+1

xk

xk+1


xk

f (x) dx =

a

a

f (x) dx + ··· +

f (x) dx + ··· +

b

f (x) dx =

b − a
2n



n−1

k=1

f (a) + f (b) + 2

f (xk)

+ R(ξ),

|R(ξ)| ≤ M

(b − a)2
12n2

, M = máx
x∈[a,b]

|f(x)|.

donde

a

Problema 4 Utilizando las fórmulas (6) y (7) demostrar las fórmulas (10) y (11).

Problema 5

Prueba la fórmula (6) y (7).

Implementación con Maxima.

En este caso volvemos a definir la partición:

(%i11)

kill(all)$
a:0;b:1;x[0]:a;Nu:20;
x[n]:=x[0]+n*(b-a)/Nu;

(%o1)
(%o2)
(%o3)
(%o4)
(%o5)

0
1
0
20
x[n]:=x[0]+(n*(b−a))/Nu

y definimos la función y la suma numérica

(%i6)

define(f(x),x^2);
tra: (f(a)+f(b)+ 2*sum(f(x[k]),k,1,Nu-1))*((b-a)/(2*Nu))$
float(%);

(%o6) f(x):=x^2
(%o8) 0.33375

Finalmente, comparamos el resultado numérico con el valor exacto

(%i9) exac:float(integrate(f(x),x,a,b));

float(abs(tra-exac));
0.33333333333333
4.166666666666763*10^−4

(%o9)
(%o10)

Método de Simpson.

(10)

(11)

b

a

f (x) dx de

El método de Simpson para calcular integrales consiste en aproximar la integral

la siguiente forma

b

a

a + b



2

f (x) dx = A f (a) + B f

+ C f (b) + R(ξ),

(12)

donde A, B, C son tales que R(ξ) es igual a cero si f (x) = 1, f (x) = x y f (x) = x2, respectivamente.
Es decir si sustituimos en (12) la función f por cualquiera de las funciones f (x) = 1, f (x) = x o
f (x) = x2, la fórmula es exacta, o sea R(ξ) = 0. Esto es equivalente a aproximar el área debajo de f
por una parabola (ver figura 3). Este método es un caso particular del método de Newton-Cotes.

4

Figura 3: Aproximación de una integral por el método de Simpson. A la izquierda vemos el área bajo
la curva que queremos calcular. A la derecha, la aproximación usando el método de los trapecios que
equivale a encontrar el área bajo una parábola.

Problema 6 Sustituyendo f (x) = 1, f (x) = x y f (x) = x2 en (12) encontrar un sistema de ecuacio-
nes para las incógnitas A, B, C y demostrar entonces que (13) se puede escribir de la forma

f (x) dx =

b − a
6

f (a) +

4(b − a)

6

f

b − a
6

+

f (b) + R(ξ).

(13)

a + b



2

b

a

Si f es cuatro veces derivable y todas sus derivadas son continuas en [a, b] entonces se puede demostrar
que R(ξ) se expresa de la forma

R(ξ) =

(b − a)5
2880

f (4)(ξ),

ξ ∈ [a, b].

Problema 7 Demostrar la fórmula anterior. Para ello seguir los siguientes pasos.

1. Comprobar que la función F (x, t), con x = a+b

2 , definida por

x+t

x−t

F (x, t) =

f (ξ)dξ − t
3

[f (x − t) + 4f (x) + f (x + t)] ,

(14)

(15)

es continua y tres veces diferenciable para todo t ∈ [0, b−a
además F (x, 0) = F (x, 0) = F (x, 0) = 0.

2 ], y F (x, t) = − t

3 [f(x+t)−f(x−t)],

2. Probar que F (x, t) es tal que existen dos números reales m y M (¿quiénes son dichos números?)

tales que

y deducir de aquí que

2
3

mt2 ≤ F (x, t) ≤ 2
3

M t2,

1
90

mt5 ≤ F (x, t) ≤ 1
90

M t5.

3. Finalmente, substituyendo t = b−a

2 , deducir el resultado deseado.

b

Al igual que en los casos anteriores vamos aproximar la integral

f (x) dx con mejor exactitud

dividiendo el intervalo [a, b] en 2n puntos de la forma

a

[a, b] = [a, x1] ∪ [x1, x2] ∪ ··· ∪ [x2n−2, x2n−1] ∪ [x2n−1, b],

donde

xk = a +

b − a
2n

k,

k = 0, 1, 2, ..., 2n,

x0 = a, x2n = b.

5

yxf(x)yxf(x) b

b

x2



f (x) dx =

b − a
6n

x2k+2

n

b



n−1

Apliquemos ahora la fórmula de Simpson (13) para cada subintervalo [x2k, x2k+2], k = 0, 1, ..., n − 1,
o sea, escribamos la integral original como la suma de las integrales

f (x) dx =

f (x) dx + ··· +

f (x) dx + ··· +

a

a

x2k

f (x) dx.

x2n−2

y apliquemos el método de Simpson a cada uno de los sumandos. Nótese que los intervalos siguen
teniendo una longitud x2k+2 − x2k = b−a

n igual que antes. Esto nos conduce a la expresión

a

donde

f (a) + f (b) + 4

f (x2k−1) + 2

f (x2k)

+ R(ξ),

(16)

k=1

k=1

|R(ξ)| ≤ M

(b − a)5
2880n4 , M = máx
x∈[a,b]

|f (4)(x)|.

(17)

Problema 8 Utilizando las fórmulas (13) y (14) demostrar las fórmulas (16) y (17).

Implementación con Maxima.

En este caso la partición es de 2n puntos. Así definimos

(%i11)

kill(all)$
a:0;b:1;x[0]:a;Nu:10;
x[n]:=x[0]+n*(b-a)/(2*Nu);

(%o1) 0
(%o2) 1
(%o3) 0
(%o4) 10
(%o5) x[n]:=x[0]+(n*(b−a))/(2*Nu)

y definimos la función y la suma numérica

(%i6) define(f(x),x^2);

simp: (f(a)+f(b)+ 4*sum(f(x[2*k-1]),k,1,Nu)+

2*sum(f(x[2*k]),k,1,Nu-1))*((b-a)/(6*Nu))$

float(%);

(%o6) f(x):=x^2
(%o8) 0.33333333333333

Finalmente, comp