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

Imágen de pdf 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 CASgráfica de visualizaciones

Publicado el 13 de Julio del 2018
347 visualizaciones desde el 13 de Julio del 2018
579,1 KB
13 paginas
Creado hace 2a (26/06/2017)
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
  • Links de descarga
http://lwp-l.com/pdf12500

Comentarios de: Integración numérica y resolución de algunas ecuaciones diferenciales ordinarias con Maxima CAS (0)


No hay comentarios
 

Comentar...

Nombre
Correo (no se visualiza en la web)
Valoración
Comentarios
Es necesario revisar y aceptar las políticas de privacidad