*
“FACULTAD DE CIENCIAS MATEMATICAS”
METODOS NUMERICOS I-II PROF:LUCIO AVILIO MALASQUEZ RUIZ.
ERROR ABSOLUTO (𝜺𝒂 ) Es la diferencia entre un valor exacto y una d sus aproximaciones. Puede ser positivo o negativo, según si la medida es superior al valor real o inferior (la resta sale positiva o negativa). Tiene unidades, las mismas que las de la medida. Ahora sea: P = ln(5) Entonces: P= ln(5) = 1.609437912
Y sea sus siguientes aproximaciones: P1 = 1.60123456 P2 = 1.609413131 P3 = 1.6094372
DEFINICION DE CIFRAS SIGNIFICATIVAS EXACTAS: Se define por la siguiente relación: 𝜀𝑎 = 0.5 × 10𝑚−𝑛+1 ; Donde: m=es la descomposición polinómica. n= numero de cifras significativas exactas.
DEFINICION DE CIFRAS DECIMALES EXACTAS: Se define por la siguiente relación: 𝜀𝑎 = 0.5 × 10−𝑡 ; Donde: t=numero de cifras decimales exactas. Ahora sea: P = ln(5) Entonces: P= ln(5) = 1.609437912 Y sea sus siguientes aproximaciones: P1 = 1.60123456 P2 = 1.609413131 P3 = 1.6094372
Donde: m = 0 Ahora hallamos sus Ea de cada uno. Ea1 = | p – p₁ | = 0.008203 = 0.008 = 0.8x10¯³ = 0.001 = 0.1x10¯² Ea₂ = |p - p₁ | = 0.000025 = 0.00003 = 0.3x10¯⁴ Ea₃ = |p - p₃ | = 0.0000012 = 0.00001 = 0.1x10¯⁵
Ahora hallaremos el número de cifras significativas exactas: Para p₁: 0.1x10¯² ≤ 0.5x10¯ⁿ⁺¹ -2 = -n+1 n = 3 ---> p₁ tiene 3 cifras significativas exactas.
Ahora hallaremos el número de cifras decimales exactos: Para p₁: 0.1x10¯² ≤ 0.5x10¯™ K=2 Entonces p₁ tiene 2 cifras decimales exactas. Para p₂: 0.3x10¯⁴ ≤ 0.5x10¯™ K=4 Entonces p₂ tiene 4 cifras decimales exactas. Para p₃: 0.1x10¯⁵ ≤ 0.5x10¯™ K=5 Entonces p₃ tiene 5 cifras decimales exactas.
TEOREMA DE BOLZANO (TB) Sea la ecuación no lineal f(x) =0; Donde f(x) es función no trascendente (trigonométrica, exponencial, logarítmica o polinomial), 𝑥 ∈ [𝑎, 𝑏] Si 𝑓(𝑎)𝑓(𝑏) < 0, 𝑒𝑛𝑡𝑜𝑛𝑐𝑒𝑠 𝑒𝑥𝑖𝑠𝑡𝑒 𝑙𝑎 𝑟𝑎𝑖𝑧 𝑜 𝑠𝑜𝑙𝑢𝑐𝑖𝑜𝑛 𝑥 ∗ ∈< 𝑎, 𝑏 >/. 𝑓(𝑥 ∗ ) ≡ 0
Aplicación del Teorema de Bolzano: Sea f(x) = x.cos(x) – 1 x
f(x)
-2
-
-1
-
0
-
1
-
2
+
Como f(1).f(2) < 0 ---> existe x* que pertenece a [1,2]
LOCALIZACIÓN DEL INTERVALO DONDE EXISTE LA RAÍZ
POR MÉTODO GRAFICO
POR MÉTODO GRAFICO
•Paso 4: Proyectar hacia el eje x, los puntos de intersección de
MÉTODO DE BISECCIÓN Dado la ecuación no lineal f(x) = 0 /. Existe la raíz X* ∈ [a, b] por T.B. El método de bisección consiste en hallar el promedio simple de cada intervalo
EJEMPLO
MÉTODO DE REGULA FALSI O MÉTODO DE FALSA POSICIÓN Es un método similar al método de bisección en la que en vez de hallar el promedio simple de dos intervalos --> X*. El método de R.F. determina el promedio ponderado de los intervalos --> X*.
Caso contrario ir al P- 2
Ejemplo
MÉTODO DE PUNTO FIJO O MÉTODO DE APROXIMACIONES SUCESIVAS
EJEMPLO
MÉTODO DE LA SECANTE
EJEMPLO
MÉTODO DE NEWTON – RAPHSON El Método de Newton-Raphson es ampliamente utilizado para encontrar las raíces de la ecuación f(x)=0, ya que converge rápidamente, la contra es que uno debe conocer la derivada de f(x) y se necesita una aproximación inicial muy cercana a la raíz. Se requiere que f(x) sea doblemente continua y diferenciable en [a,b].
Ejemplo:
Sea f(x)=-3x ; donde x* ∈ [0,1] Primero analizamos su convergencia:
Ejemplo aplicativo
INTERPOLACIÓN Supongamos que se conoce f0 , f1, f2, …….fn valores correspondientes a X0, X1, X2, ….., Xn valores independientes de una variable independiente X.( X0
a) Interpolación directa.- consiste en que dado un valor XP diferente de los Xi pero correspondido entre X0 y Xn, se desea hallar el valor de su imagen fP
a) Interpolación inversa.- consiste en que dado el valor de la imagen fP se desea hallar el valor XP que genera dicha imagen.
1.- INTERPOLACIÓN DIRECTA LINEAL XP = valor a interpretar
f
FP = f (XP) imagen de XP F(X)
(X1,f1) (X0,f0) (XP,fP)
Xo
XP n
X1
X
INTERPOLACIÓN DIRECTA DE NEWTON -PROGRESIVO Y NEWTON – REGRESIVO. Para un conjunto de (n+1) puntos igualmente espaciados Interpolación consiste, en dado un valor no considerado 𝑥𝑝 , en la tabla, se debe hallar su imagen 𝑓𝑝 . Hay dos tipos de Interpolación. A) Interpolación Directa. Consiste en que dado 𝑥𝑝 se debe hallar su imagen 𝑓𝑝 . B) Interpolación Inversa. Consiste en que dado el valor de la imagen 𝑓𝑝 , se debe hallar el valor 𝑥𝑝 .
1) Interpolación Directa de Newton Progresivo(IDNP) Se utiliza cuando se desea interpolar un valor xp dado al principio de la tabla o 1er sector de la tabla. Utiliza la siguiente fórmula.
𝑓𝑝 = 𝑓0 + P∆𝑓0 + P P−1)…(P−(n−1) n!
P(P−1) 2!
∆2 𝑓0 +
P P−1 (P−2) 3!
∆3 𝑓0 + ⋯ +
∆n 𝑓0
Donde P=xp -x0 h , xp∈x0, x1 , P∈ <0,1> Utiliza la siguiente Tabla de omisión de Términos (T.0.T) ∆2 Newton Progresivo Newton Regresivo
4
∆3 8
∆4 12
∆5 16
Observación: Si |∆2 𝑓𝑖 |<4 Se tiene Inter. Directa Lineal y se utiliza la fórmula de N.P. hasta la 1𝑟𝑎 Diferencia, o sea; 𝑓𝑝 = 𝑓0 +𝑃∆𝑓0 . Si |∆3 𝑓𝑖 |<8 Se tiene Inter. Directa No Lineal (IDNL) y se utiliza la fórmula de N.P. hasta la 2𝑑𝑎 Diferencia, o sea; 𝑓𝑝 = 𝑓0 +𝑃∆𝑓0 + P(P−1) 2 ∆ 𝑓0 2!
.
Si |∆4 𝑓𝑖 |<12 Se utiliza IDNL y en la formula de NP se utiliza hasta la 3𝑟𝑎 diferencia, osea 𝑓𝑝 = 𝑓0 +𝑃∆𝑓0 +
P(P−1) 2 P P−1 (P−2) 3 ∆ 𝑓0 + ∆ 𝑓0 . 2! 3!
(4) Tanto en (1),(2) y (3) solo se considera las cifras significativas.
*
2) Interpolación de Newton Regresivo (NR) Se utiliza cuando se quiere interpolar un valor en la parte final de la tabla. Su fórmula es
𝑓𝑝 = 𝑓0 + P𝛻𝑓0 +
P(P+1)
𝛻 2 𝑓0 +
2! P P+1)…(P+(n−1) n!
P P+1 (P+2) 3!
𝛻 3 𝑓0 + ⋯ +
𝛻 n 𝑓0
Donde 𝑃 ∈< −1,0 > 𝑦 𝑥𝑝 = ,𝑥−1 , 𝑥0 -
Interpolación Directa Central
Interpolación de Stirling
Interpolación de Bessel
Interpolación de Everett
Ejemplo
INTERPOLACIÓN PARA EL CASO DE DATOS NO EQUIDISTANTES
POLINOMIO DE APROXIMACIÓN DE LAGRANGE
Sean
(X0, f0), (X1, f1),…, (Xn,fn)
(n+1) ptos
Entonces existe un polinomio de grado ≤ 𝑛 que para por dichos puntos
(X0
OBS: en el primer termino falta (X-X0) , en el 2º termino (X-X1), y asi sucesivamente en el ultimo termino falta (X-Xn) esto es una cualidad de dicho polinomio. Como f(x) debe contener a los puntos dados
Si X = X0
; f(X0)= a0(X0-X1)(X0-X2)……..(X0-Xn)
a0= Si X = X1
f(o) X0−X1 X0−X2 X0−X3 ………(X0−Xn)
; f(X1)= a1(X1-X0)(X1-X2)……..(X1-Xn)
a1=
f(1) X1−X0 X1−X2 X1−X3 ………(X1−Xn)
. .
. Si X = Xn
; f(Xn)= an(Xn-X0)(Xn-X1)(Xn-X2)……..(Xn-Xn-1)
an=
f(n) Xn−X1 Xn−X2 Xn−X3 ………(Xn−Xn−1)
TEOREMA.- Sean
(X0, f0), (X1, f1),…, (Xn,fn) para los puntos (n+1) y además (X0
INTERPOLACIÓN Y APROXIMACIÓN DE LAGRANGE
Polinomio de Lagrange
Resolver un polinomio de Lagrange
Ejemplo
METODOS DIRECTOS DE SISTEMAS DE ECUACIONES LIENALES
FACTORIZACIÓN MATRICIAL “LU” (Método de Crout – Doolitle) Consiste en factorizar una matriz cuadrada “A” en un producto “LU”. Esto es: 𝐴 = 𝐿. 𝑈
Donde: A: es la matriz a factorizar L : es una matriz triangular inferior, cuyos elementos de la diagonal principal son iguales a 1. Se llama “L” porque viene de la palabra inglesa “low”, que significa “bajo”. U : es una matriz triangular superior, cuyos elementos se hallan por el método de la eliminación gaussiana. Se llama “U” porque viene de la palabra inglesa “up”, que significa “arriba”.
PRIMER METODO DIRECTO:
MÉTODO CROUT-DOOLITLE.
Ejemplo
SEGUNDO METODO DIRECTO METODO DE CHOLESKY Resolver por método de Cholesky
Ejemplo
MÉTODO TRIDIAGONAL PARA SISTEMAS DE ECUACIONES LINEALES.
SISTEMA DE ECUACIÓN
ALGORITMO TRIDIAGONAL
EJEMPLO
MÉTODO PENTADIAGONAL PARA SISTEMAS DE ECUACIONES LINEALES
⋯ ⋱ ⋯
⋮ 0
0 ⋮
A
𝑥1 𝑥2 ⋮ ⋮ 𝑥𝑛
=
x
𝑥1 𝑥2 ⋮ ⋮ 𝑥𝑛
Ec (1) Ec (2) ⋮ Ec (n)
b
Algoritmo Pentadiagonal
P-1 Dado el sistema Ax=b , expresarlo en la forma A𝑥 =b y definir 𝑥1 =C ∧ 𝑥2 =D
Donde C ∧ D son constantes arbitrarias.
P-2 De Ec (1) despejar 𝑥3
De Ec (2) despejar 𝑥4 ⋮
De Ec (n-2) despejar 𝑥𝑛 De Ec (n-1) despejar 𝑥𝑛+1 De Ec (n) despejar 𝑥𝑛+2 Como no existe 𝑥𝑛+1 se hace lo siguiente : 𝑎𝑛−1,𝑛−2 𝑥𝑛−2 + 𝑎𝑛−1,𝑛−1 𝑥𝑛−1 = 𝑏𝑛−1 + 𝑅1
Como no existe 𝑥𝑛+2 se hace lo siguiente : 𝑎𝑛,𝑛−1 𝑥𝑛−1 + 𝑎𝑛,𝑛 𝑥𝑛 = 𝑏𝑛 + 𝑅2
Donde 𝑅 =
𝑅1 es un vector residual. 𝑅2
Si 𝑅 = 0 → 𝑥1 = (𝑥1 , 𝑥2 , … , 𝑥𝑛 )𝑡 es solución del sistema A.x = b Si 𝑅 ≠ 0 ⟶ 𝑥1 𝑒𝑠 𝑠𝑜𝑙𝑢𝑐𝑖ó𝑛 𝑑𝑒𝑙 𝑠𝑖𝑠𝑡𝑒𝑚𝑎 𝐴. 𝑥1 = 𝑏 + 𝑅
P-3 : La primera solución homogénea del sistema A.x = b se debe expresar como
𝐴𝑥 = 𝜃, 𝑑𝑜𝑛𝑑𝑒 𝜃 𝑒𝑠 𝑒𝑙 𝑣𝑒𝑐𝑡𝑜𝑟 𝑛𝑢𝑙𝑜; luego definir 𝑥1 = 𝐶 𝑦 𝑥2 = 𝐷 en la cual se considera C y D como constantes arbitrarias, finalmente se procede similar al paso anterior (P-2) llegando a la siguiente solución 𝑎𝑛−1,𝑛−2 𝑥𝑛−2 + 𝑎𝑛−1,𝑛−1 𝑥𝑛−1 = 𝜃 + 𝑠1 𝑠1 𝑎𝑛,𝑛−1 𝑥𝑛−1 + 𝑎𝑛,𝑛 𝑥𝑛 = 𝜃 + 𝑠2 en donde 𝑆 = .𝑠 / es un vector residual. 2
Si 𝑆 = 0 → 𝑥2 = (𝑥1 , 𝑥2 , … , 𝑥𝑛 )𝑡 es solución del sistema A.x = b Si 𝑆 ≠ 0 ⟶ 𝑥2 𝑒𝑠 𝑠𝑜𝑙𝑢𝑐𝑖ó𝑛 𝑑𝑒𝑙 𝑠𝑖𝑠𝑡𝑒𝑚𝑎 𝐴. 𝑥2 = 𝜃 + 𝑆
P-4 : La segunda solución homogénea del sistema A.x = b se debe expresar como 𝐴𝑥 = 𝜃, 𝑑𝑜𝑛𝑑𝑒 𝜃 𝑒𝑠 𝑒𝑙 𝑣𝑒𝑐𝑡𝑜𝑟 𝑛𝑢𝑙𝑜; luego definir 𝑥1 = 𝐶 𝑦 𝑥2 = 𝐷 en la cual se considera C y D como constantes arbitrarias, finalmente se procede similar al paso P-2 llegando a la siguiente solución.
𝑎𝑛 −1,𝑛 −2 𝑥𝑛 −2 + 𝑎𝑛 −1,𝑛 −1 𝑥𝑛 −1 = 𝜃 + 𝑇1 𝑎𝑛 ,𝑛 −1 𝑥𝑛 −1 + 𝑎𝑛 ,𝑛 𝑥𝑛 = 𝜃 + 𝑇2 en donde 𝑇 =
Si 𝑇 = 0 → 𝑥3 = 𝑥1 , 𝑥2 , … , 𝑥𝑛
𝑡
𝑇1 𝑇2
es un vector residual.
es solución del sistema A.x = b
Si 𝑇 ≠ 0 ⟶ 𝑥3 𝑒𝑠 𝑠𝑜𝑙𝑢𝑐𝑖ó𝑛 𝑑𝑒𝑙 𝑠𝑖𝑠𝑡𝑒𝑚𝑎 𝐴. 𝑥3 = 𝜃 + 𝑇 P-5
: Finalmente se llegará al siguiente sistema 𝐴. 𝑥1 = 𝑏 + 𝑅 ………(1) 𝐴. 𝑥2 = 𝜃 + 𝑆
………(2)
𝐴. 𝑥3 = 𝜃 + 𝑇 ……….(3)
Luego hacemos (1) – 𝛼(2) − 𝛽(3) y se llega a lo siguiente
A (𝑥1 − 𝛼(2) − 𝛽(3)) = 𝑏 + (𝑅 − 𝛼𝑆 − 𝛽𝑇 )
X
𝜃
Se busca luego una solución en donde R −𝛼𝑆 − 𝛽𝑇 = 𝜃 / 𝑋 = 𝑥1 − 𝛼𝑥2 − 𝛽𝑥3 Entonces la solución del sistema A.x = b estará dado por
=
1
−
2
−
3
Ejemplo Resolver el siguiente sistema pentadiagonal
2𝑋1+ 3𝑋2+ 𝑋3
= 8 Ec(1)
3𝑋1+ 2𝑋2 + 4𝑋3 + 𝑋4
= 15 Ec(2)
𝑋1+ 4𝑋2 + 𝑋3 + 4𝑋4 + 2𝑋5
= 13 Ec(2)
𝑋2 + 4𝑋3 + 2𝑋4 + 𝑋5
= 19 EC(4)
2𝑋3 + 𝑋4 + 7 𝑋5
= 15 EC(5)
PASO DEL ALGORITMO
P-1: Expresar el sistema como A. 𝑋= b 2𝑋1 + 3𝑋2 + 𝑋3
= 8 Ec(1)
3𝑋1 + 2𝑋2 + 4𝑋3 + 𝑋4
= 15 Ec(2)
𝑋1 + 4𝑋2 + 𝑋3 + 4𝑋4 + 2𝑋5
= 13 Ec(3)
𝑋2 + 4𝑋3 + 2𝑋4 + 𝑋5
= 19 Ec(4)
2𝑋3 + 𝑋4 + 7𝑋5
= 15 Ec(5)
Sea 𝑋1 = 0
⋀ 𝑋2 = 1
cte arbitrario
P-2: DE Ec(1) despejar 𝑋3 → 𝑋3 = 5 pues 0 + 3(1) + 𝑋3 = 8 DE Ec(2) despejar 𝑋4 → 𝑋4 = −7 pues 0 + 2(1) +4(8) + 𝑋4 = 15 DE Ec(3) despejar 𝑋5 → 𝑋5 = 16 pues 0 + 4(1) + 5 + 4(−7) + 2𝑋5 = 13
DE Ec(4) despejar 𝑋6 ; como ∄ 𝑋6 , hacemos lo sgte. : 𝑋2 + 4𝑋3 + 2𝑋4 + 𝑋5 = 19 + 𝑅1
1 + 4(5) + 2(−7) + 16 = 19 + 𝑅1 → 𝑅1 = −45 De la Ec(5) despejar 𝑋7 ; como 𝑋7 ∄ hacemos 2𝑋3 + 𝑋4 + 7𝑋5 = 15 + 𝑅2 2(5) + (−7) +7(16) = 15 + 𝑅2 → 𝑅2 = 100
R=
−45 100
, como R ≠ 0; A. 𝑥1 = b + R
Se tiene:
𝑥1 =
𝑋1 𝑋2 𝑋3 𝑋4 𝑋5
=
0 1 5 −7 16
ó 𝑋1 = 𝑋1 , 𝑋2 , 𝑋3 , 𝑋4 , 𝑋5
𝑡
P-3: Primera solución homogénea del sistema A.x = b, expresarlo Como A.𝑥= 𝜃 2𝑋1 + 3𝑋2 + 𝑋3
= 0 Ec(1)
3𝑋1 + 2𝑋2 + 4𝑋3 + 𝑋4
= 0 Ec(2)
𝑋1 + 4𝑋2 + 𝑋3 + 4𝑋4 + 2𝑋5
= 0 Ec(3)
𝑋2 + 4𝑋3 + 2𝑋4 + 𝑋5
= 0 Ec(4)
2𝑋3 + 𝑋4 + 7𝑋5
= 0 Ec(5)
Sea 𝑋1 = 10 ⋀ 𝑋2 = 20 De Ec(1) despejar 𝑋3 ; → 𝑋3 = −80
pues 2(10) + 3(29) + 𝑋3 = 0
De Ec(2) despejar 𝑋4 ; → 𝑋4 = 250
𝑥2 =
pues 3(10) +2(20) + 4(−8) + 𝑋4 = 0 De Ec(3) despejar 𝑋5 ; → 𝑋5 = −505 pues 10 + 4(20) + (−80) + 4(250) + 2𝑋5 =0 De Ec(3) despejar 𝑋6 ; como ∄ 𝑋6 hacemos 𝑋2 + 4𝑋3 + 2𝑋4 + 𝑋5 = 0 + 𝑆1
20 + 4(−80) + 2(250) + (−505) = 0 + 𝑆1 ; → 𝑆1 = 1445 De Ec(5) despejar 𝑋7 ; como ∄ 𝑋7 hacemos 2𝑋3 + 𝑋4 + 7𝑋5 = 0 + 𝑆2 ; → 𝑆2 = −3445 S=
1445 −3445 2(−80)+ (250) + 7(−505)= 0 + 𝑆2 = 0
10 20 −80 250 −505
P-4: Segunda solución homogénea del sistema A.x = b expresar A.𝑥= 0 2𝑋1+ 3𝑋2+ 𝑋3
= 0 Ec(1)
3𝑋1+ 2𝑋2+ 4𝑋3+ 𝑋4
= 0 Ec(2)
𝑋1 + 4𝑋2 + 𝑋3 +4 𝑋4 + 2𝑋5 = 0 Ec(3) 𝑋2 + 4𝑋3 +2𝑋4 + 𝑋5 = 0 Ec(4)
2𝑋3+ 𝑋4+ 7𝑋5 = 0 Ec(5) Sea 𝑋1 = 20 ⋀ 𝑋2 = 10 De Ec(1) despejar 𝑋3; → 𝑋3= −70
pues 2(20) + 3(10) + 𝑋3 = 0 De Ec(2) despejar 𝑋4 → 𝑋4 = 200 3𝑋1 + 2𝑋2 + 4𝑋3 + 𝑋4 = 0 3(20) + 2(10) + 4(−70) + 𝑋4 = 0 De Ec(3) despejar 𝑋5 → 𝑋5 = 395 𝑋1 + 4𝑋2 +𝑋3 + 4𝑋4 + 2𝑋5 = 0
20 + 4(10) + (−70)+ 4(200) + 2𝑋5 = 0
De Ec(4) despejar 𝑋6 como ∄ 𝑋6 , hacemos 𝑋2 + 4𝑋3 +2𝑋4 + 𝑋5 = 0 + 𝑇1 → 𝑇1 = 1925
10 + 4(−70) +2(200) + 395 = 0 + 𝑇1
De Ec(5) despejar 𝑋7 ; como ∄ 𝑋7 𝑎𝑐𝑒𝑚𝑜𝑠 2𝑋3 + 𝑋4 +7𝑋5 = 0 + 𝑇2 𝑋1
T = .𝑇𝑇1 /= 2
1995 −2705
𝑋2
; 𝑋3 = 𝑋3 𝑋4 𝑋5
A.𝑥3 = 𝜃 + T
P-5: Y se llega a lo siguiente A.𝑥1 = b + R A.𝑥2 = 0 + S A.𝑥3 = 0 + T
20 10 = −70 200 395
−45 / 100
R – 𝛼𝑆 – 𝛽𝑇 = 0
R= .
1445 / −3445
𝛼𝑆 + 𝛽𝑇 = R
S=.
𝟏𝟒𝟒𝟓 𝟏𝟒𝟒𝟓 −𝟒𝟓 /+ 𝜷 . /= . / −𝟑𝟒𝟒𝟓 −𝟑𝟒𝟒𝟓 𝟏𝟎𝟎
𝜶.
1445 𝛼+ 1925 𝛽 = −45
;
−3445𝛼 −2705 𝛽 = 100
;
𝟏𝟗𝟐𝟓 / −𝟐𝟕𝟎𝟓
T=.
𝛼 = −0.025992507 𝛽 = −0.003865364
X = 𝑋1 − 𝛼𝑋2 − 𝛽𝑋3
X=
𝑋1 𝑋2 𝑋3 𝑋4 𝑋5
=
0 1 5 −7 16
+ 0.025992507
10 20 −80 + 0.003865364 250 −505
20 10 −70 200 395
𝑋1 0.33723235 𝑋2 1.55850378 𝑋 2.65002396 X= 3 = 𝑋4 0.27119955 1.34696518 𝑋5
Comprobación :
Ec(1)
2𝑋1 + 3𝑋2 + 𝑋3 = 8
0.6744647 + 4.67551134 + 2.65002396 = 8 8
=
8 Cumple!!!
Y tambien cumple todas las ecuaciones…
NORMA DE UNA MATRIZ La Norma de una matriz 𝐴𝑛×𝑛 es un número real tal que satisface las siguientes condiciones
() ( ) 1 ( ) ( )
>0 =
= 0,
=0
, +
≤
−
≤
=
,
−1 =
+
Principales Normas
-Norma
“m”
-Norma “l”, -Norma “k”,
o
→
Norma = =
2 ,
=
Ejemplo Sea 𝐴 =
−1 −2 3 4 5 6 7 8 9
𝐴
𝑚
= max* −1 + −2 + 3,4 + 5 + 6,7 + 8 + 9+ = max*6,15,24+ = 24
𝐴
𝑙
= max* −1 + 4 + 7, −2 + 5 + 8,3 + 6 + 9+ = max*12,15,18+ = 18
𝐴
𝑘
=
(−1)2 + (−2)2 + 32 +. . +82 + 92 =
𝑥1 𝑥2 Para el vector X = : . 𝑥𝑛
→
𝑋
𝑚
285 = 16.9
= max*𝑥 +
𝑋
𝑙
𝑋
𝑘
= 𝑥1 + 𝑥2 +. . +|𝑥𝑛 | =
|𝑥1 |2 + |𝑥2 |2 + |𝑥3 |2
SOLUCION ITERATIVA DE SISTEMAS DE ECUACIONES LINEALES
PRIMER METODO: METODO DE JACOBI Dado el sistema AX=b………… (1) Despejamos X de la ecuación 1 obteniendo un sistema equivalente de la forma: X = β + αX……… (2) De la siguiente manera.
Ec.(1) a11X1+a12X2+ …………………………………………………….. +a1mXm = b1 Ec.(1) a21X1+a22X2+ …………………………………………………….. +a2mXm = b2 Ec.(1) a31X1+a32X2+ …………………………………………………….. +a3mXm = b3 . . . Ec.(1) am1X1+am2X2+ …………………………………………………….. +ammXm = b1
11
Donde:
A=
⋮ 1
⋯ ⋱ ⋯
1
1
⋮
b = ⋮ y si aii ≠ 0
1
X= ⋮
Despejamos X; de la Ec. (i);
i=1,2,3, …. ,m
X = β + αi-1X
X1 = X2 =
𝑏1 𝑎 11 𝑏2 𝑎 22
−
𝑎 12
−
𝑎 21
𝑎 11 𝑎 22
𝑋2 −
𝑎 13
𝑋1 −
𝑎 23
. . .
𝑎 11 𝑎 22
𝑋3 − … … … … −
𝑎1m
𝑋3 − … … … … −
𝑎2m
𝑎11 𝑎22
𝑋m 𝑋m
………. (2)
Donde:
β = (β1,β2, … … … ,βm)t
β = bi/aii ; aii ≠ 0 ;
α = αij = -aij/aii
0 21
α=
⋮
12 0
1
2
⋯ ⋱ ⋯
1 2 ⋮ 0
El sistema (2) anterior se puede expresar en la siguiente forma X = β + αX
X
=
β
+
αX
0 𝛼12 𝑋1 𝑋1 ⋮ = ⋮ + 𝛼21 ⋮ 0 𝑋m 𝑋m 𝛼m1 𝛼m2
⋯ ⋱ ⋯
𝛼1m 𝑋1 𝛼2m ∗ ⋮ ⋮ 𝑋m 0
… (2)
El Sistema (2) sugiere Jacobi la siguiente relación de recurrencia
X(k+1) = β + αX(k) ,
k=0,1,2, … … …
Ó
…………………….. (3)
X(k) = β + αX(k-1),
k=1,2, … … …
De la relación (3) se obtiene la sucesión {Xk}∞k=0 tomando como vaor inicial X(0) arbitrario, que generalmente X(0)=0 ó X(0)=β ó β=1 Obs.
X(k+1) = (X1(k+1),X2(k+1), … … … … … … … … … … … … … , Xm(k+1))t X(0) = (X1(0),X2(0), … … … … … … … … … … … … … … , Xm(0))t
ALGORITMO DE JACOBI P-1
Dado el Sistema Ax = b …………………………(1) Expresarlo en el sistema equivalente X = β + αx ……………………(2)
P-2 Tomando como solución inicial X(0) arbitrario generar la sucesión {X(k)} → X(*) mediante la relación de recurrencia: X(k+1) = β + αx(k)
,
k=0,1,2, … …. …. …
,
k=1,2, … … …
Ó X(k) = β + αX(k-1) P-3
Dejar de iterar si
(|| X(k) – X(k-1) || )/ || X(k) || <= ε = 10-m C.C. ir al P-2
NOTACIÓN MATRICIAL DEL MÉTODO DE JACOBI Sea el sistema Ax = b Donde:
A= 11
12
21
22
⋮
⋮
1
2
… … ⋱ …
1 2
⋮
La matriz A se le puede descomponer en la forma A = D +L + U, donde
0
11
=
0 0
22
⋮ 0
Matriz Diagonal
⋯ ⋱ ⋯
0 0 ⋮
0 =
21
0 0
⋯
1
…
(
⋮
Matriz Triangular inferior
⋱ 1)
0 0 ⋮ 0
0 = 0 0
12
⋮
0 0
⋯
1
⋱ ⋯
⋮ 0
Matriz Triangular superior
2
Así el sistema Ax = b se le puede expresar como: (D + L + U)X = b DX + (L + U) X = b DX = b – (L + U) X X = D-1b – D-1(L + U) X X = D-1b + [-D-1(L + U)] X
→
β = D-1b
Si el método de Jacobi es
X = β + αX
^
α = -D-1(L + U)
Matricialmente es: → X = D-1b – D-1(L + U) X Su relación de recurrencia es
X (k+1)= D-1b – D-1 (L + U) X (k)
,
k=0, 1, 2 …
EJEMPLO Sea el siguiente sistema Ec (1) 20x1 + 5x3 =2 Ec (2) x1 + 20x2 + 2x3 = 4 Ec (3) x1 + 9x2 + 20x3 = 6
por Jacobi verificar su convergencia
CONVERGENCIA DE JACOBI {x (x) } X*
Si ||α|| < = 1
…………….. (i)
Observación Para que se cumpla (i) Es necesario que A del sistema original Ax = b sea diagonalmente dominante, es decir aii >= ∑ |aij| de su fila y de su columna Así: a11 = 20 > = 0 + 5 de su fila a11 = 20 > = 1 +1 de su columna a22 = 20 > = 1 + 2 de su fila a22 = 20 > = 0 + 9 de su columna Igual para a33
x1 = x2 = x3 =
2 20
+0+0-
4
-
20 6
-
20
1 20 1 20
5 20
𝑥3 2
𝑥1 + 0 𝑥1 -
9 20
20
𝑥3
𝑥2 + 0
… …. . …………………. x = β+
αx
0 α=
1 − 20 1 − 20
5
0 − 20 0
− 9
2 20
− 20 0
0 0 − 0.2 = −0.05 0 − 0.1 −0.05 − 0.45 0 5
1
1
1
9
|| α ||∞ = máx{ 0 + 0 + | − 20 | , | − 20 | + 0 + | − 20 | , | − 20 | + | − 20 | + 0 } || α ||∞ = máx{0.25, 0.15, 0.5 } || α ||∞ = 0.5 < = 1 {
} x* por jacobi
por jacobi obtener una solución con € = Si la relación de jacobi es
=β+α
, k = 0, 1,2…
Para k = 0 Interacción inicial =β+α
Observación
Sea
es arbitraria
=
=
=
+
….. β =β=
……………………. … +
α
=
Para k = 1 Primera iteración
=
=
=β+α
+
Donde = 0.1 +
= 0.1 + (0)(0.1) + (0)(0.2) + (-0.25)(0.3) = 0.025\ = 0.2 +
= 0.2 + (-0.05)(0.1) + (0)(0.2) + (-0.1)(0.3) = 0.165 = 0.205
Para k = 2 Segunda iteración
=
=
=β+α
+
….. β
=
………………….... ……… α
…………
k =3 =β+α
Tercera iteración
=
=
+
….. β
=
………………….... ……… α
Verificamos si se llego a la solución
=
= = = x* con € <
= 0.03289625 < =
…………
SEGUNDO METODO: METODO DE GAUSS- SEIDEL También determina la solución del sistema 𝐴𝑥 = 𝑏 iterativamente. De la relación matricial del sistema 𝐴𝑥 = 𝑏 : (𝐷 + 𝐿 + 𝑈)𝑥 = 𝑏 → 𝐷𝑥 + 𝐿𝑥 + 𝑈𝑥 = 𝑏 𝐷𝑥 = 𝑏 − 𝐿𝑥 − 𝑈𝑥 → 𝑥 = 𝐷 −1 𝑏 + (−𝐷−1 𝐿)𝑥 + (−𝐷−1 𝑈)𝑥 … (𝜃) De (𝜃) se obtiene la relación matricial de G-S , siguiente: ( +1)
=
−1
+ (−
−1
)
( +1)
+ (−
−1
)
( )
Observación: *Si 𝛼 = −𝐷−1 (𝐿 + 𝑈) →
−𝐷−1 𝐿 𝑒𝑠 𝑙𝑎 𝑀𝑎𝑡𝑟𝑖𝑧 𝑇𝑟𝑖𝑎𝑛𝑔𝑢𝑙𝑎𝑟 𝐼𝑛𝑓𝑒𝑟𝑖𝑜𝑟 𝑑𝑒 𝛼 −𝐷−1 𝑈 𝑒𝑠 𝑙𝑎 𝑀𝑎𝑡𝑟𝑖𝑧 𝑇𝑟𝑖𝑎𝑛𝑔𝑢𝑙𝑎𝑟 𝑆𝑢𝑝𝑒𝑟𝑖𝑟𝑜 𝑑𝑒 𝛼
*La relación (𝜃) se puede obtener al igual que Jacobi del sistema 𝐴𝑥 = 𝑏 , de la ecuación 𝑖despejar la variable 𝑥𝑖 , para obtener la Matriz 𝛼.
ALGORITMO DEL MÉTODO DE GAUSS-SEIDEL:
Paso1: Dado el sistema 𝐴𝑥 = 𝑏 obtener su sistema 𝑋 = 𝛽 + 𝛼𝑥. Paso 2: Para un punto inicial arbitrario 𝑥 (0) generar la sucesión 𝑥 (𝑘) → 𝑥 (∗) mediante la siguiente relación: 𝑥 (𝑘+1) = 𝐷−1 𝑏 + (−𝐷−1 𝐿)𝑥 (𝑘+1) + (−𝐷−1 𝑈)𝑥 (𝑘) ; 𝐾 = 0, 1, 2, … Paso 3: Dejar de iterar si
𝑥 (𝑘 ) −𝑥 (𝑘−1 ) 𝑥 (𝑘)
≤ 𝜀 = 10−𝑛 ; caso contrario ir al paso 2.
Observación: En la convergencia del método de Gauss-Seidel también se cumple que: 𝛼 <1 →
𝑥 8𝑘) → 𝑥 (∗)
EJEMPLO 1) Dados:
𝐴=
10 2 1
3 −10 3
1 14 0 3 , 𝑏 = −5 , 𝑥 (0) = 0 10 14 0
Resuelva el sistema Ax = b por el método de Gauss-Seidel.
Solución: Utilizando Gauss-Seidel:
𝑥 (𝑘 +1) = 𝐷 −1 𝑏 + (−𝐷−1 𝐿)𝑥 (𝑘 ) + (−𝐷−1 𝑈)𝑥 (𝑘 ) Operando obtenemos la secuencia: 𝑥
(1)
7/5 39/50 = 513/500
𝑥
(3)
0.995104 = 0.995276 1.00191
𝑥
(5)
0.999792 = 0.999848 1.00007
1.0634 1.02048 0.987516
,
𝑥 (2) =
,
𝑥
(4)
1.00123 = 1.00082 0.999632
𝑥
(0)
1.00004 = 1.00003 0.999988
,
Claramente converge a la solución exacta (1, 1, 1)𝑇 . La tasa de convergencia del método de Gauss-Seidel viene dada por la norma de: 0 𝑆 = (𝐿 + 𝐷)−1 𝑈 = 0 0 Cuyas normas son: 𝐽
1=
227/500 = 0.454 y 𝐽
3/10 3/50 −6/125 ∝
1/10 −7/25 37/500
= 2/5 = 0.4.
2) Considere el siguiente sistema de ecuaciones: 3 2 1 1 𝐴 = 2 3 1 ,𝑏 = 2 1 1 3 3 ¿Puede resolver este sistema por el método de Gauss-Seidel? ¿Por qué? Si lo puede hacer, haga solo dos iteraciones a partir de la solución nula y determine la tasa numérica de convergencia. Además calcula la tasa exacta de convergencia. ¿Cuántas iteraciones necesitará para alcanzar un error absoluto de 10−5 .
Solución: El método de Gauss-Seidel es aplicable porque por que la matriz es simétrica definida positiva. Dos iteraciones conducen a: 𝑥 (0) = (0 0 0)𝑇 𝑥 (1) = (1/3 4/9 20/27)𝑇 𝑥 (2) = (−17/81 136/246 644/729)𝑇 Y la tasa de convergencia numérica la podemos calcular como (en norma infinito) 𝑥 (2) − 𝑥 46 = = 0.24 189 𝑥 (1) − 𝑥 Que se parece poco a la tasa de convergencia exacta: 𝜌((𝐿 + 𝐷)−1 𝑈) = 0.413 NOTA: Calculando con más iteraciones nos acercamos a la tasa teórica, por ejemplo: 𝑥 (10) − 𝑥 = 0.413 𝑥 (9) − 𝑥 Para alcanzar (en norma infinito) un error absoluto menor que 10−5 se requieren 13 iteraciones.
DIFERENCIA NUMERICA
Diferenciación numérica Dado una tabla de (n+1) puntos igualmente espaciados de la forma: (X0, f0), (X1, f1),… (Xn, fn) continua y diferenciable en [X0, Xn]. El problema de la Diferenciación Numérica es que dado un valor Xp ∈ [X0, Xn] se desea hallar el valor 𝑓𝑝′ tal que 𝑓 ′ (𝑋𝑝 ) = 𝑓𝑝′ , donde: P=
X p −X 0 h
… (i)
X p = X 0 + Ph f X p = fp − f(X 0 + Ph) dfp dfp dp dfp 1 f ′ Xp = = ∗ = ∗ dX p dp dX p dp h
→ f p′=
1 h
∗
df p dp
… . (ii)
Fórmula de Newton Progresivo (N P):
Si
Xp = X0 en (i) P = 0
→fórmula de NP es: fp = f0 + P∆f0 +
P(P−1)∆2 f 0 2!
+
P (P−1)(P−2)∆3 f 0 3!
+
P (P−1)(P−2)(P−3)∆4 f 0 +.. 4!
Según (ii) fp′
=
1 0∆f0 h
+
(P 2 −P)′ ∆2 f 0 2!
+
(P 3 −3P 2 +2P)′ ∆3 f 0 3!
+
(P 4 −6P 3 +11P 2 +6P)′ ∆4 f 0 +. . 1 4!
fp′
=
1 0∆f0 h
+
(2P−1)∆2 f 0 2!
+
3P 2 −6P+2 ∆3 f 0 3!
+
Si XP=X0 en (i) → P=0 en (iii) ′ fX0
=
f0′
=
1 [∆f0 h
∆2 f 0 − 2
∆3 f 0 + 3
∆4 f 0 − 4
+ ⋯]
4P 3 −18P 2 −22P−6 ∆4 f 0 4!
+1 (iii)
DIFERENCIACION DE ORDEN SUPERIOR Se desea hallar 𝑓 ′′ 𝑋𝑝 , 𝑓 ′′′ 𝑋𝑝 , en forma aproximada f ′′ X p = fp′′ además P = f
′′
fp′′ fp′′
Xp =
d(f ′ (X p ))
=
1 1 d 2 fp h h dp
=
1 d 2 fp h2 dp
dX p
=
df ′p dp
X p −X 0 h d
1 df p hd p
∗ dxp , si se sabe fp′′ = h p
En general se tiene f n p =
1 d n fp h n dp
𝒇′′𝒑 Para Newton Progresivo - hacia adelante fp′′ =
1 0∆2 f0 h2
6P 2 −18P+11 / ∆ 4 f0 … 1 12
+ (P − 1)∆3 f0 + .
Si XP=X0 → P=0 fp′′ =
1 0∆2 f0 h2
11 12
+ ∆ 3 f0 + . / ∆ 4 f0 … 1
Formula de Newton Regresivo (NR). Su formula de interpolación es: P(P + 1)∇2 f0 P(P + 1)(P + 2)∇3 f0 P(P + 1)(P + 2)(P + 3)∇4 f0 fp = f0 + P∇f0 + + + +.. 2! 3! 4! Según (ii) fp′
(P2 + P)′ ∇2 f0 (P3 − 3P2 + 2P)′ ∇3 f0 (P4 − 6P3 + 11P2 − 6P)′ ∇4 f0 1 = [∇f0 + + + +. .. h 2! 3! 4!
1 (2P + 1)∇2 f0 (3P2+ + 6P + 2)∇3 f0 (4P4 + 18P2 + 22P + 6)∇4 f0 = [∇f0 + + + +⋯ h 2! 3! 4! Si XP=0 → P=0 ′ fX0
=
f0′
=
1 0∇f0 h
∇2 f 0 + 2
∇3 f 0 + 3
∇4 f 0 + 4
+ ⋯1
Nota: Para determinar la 1ra derivada se debe tomar una diferencia más que la que indica la T O T en el caso de interpolación.
DIFERENCIACION DE ORDEN SUPERIOR 𝐟 ′′ PARA NEWTON REGRESIVO
6P 2 +18P +11 / ∇ 4 f0 … 1 12
1
fP′′ = h 2 0∇2 f0 + (P + 1)∇3 f0 + . 1
11
P = 0 → f0′′ = h 2 0∇2 f0 + ∇3 f0 + .12 / ∇4 f0 … 1 Ejemplo: En la siguiente tabla determinar (a)𝑓0′ y (b)𝑓 1 (0.055) por NP X 0.00 0.05 0.10 0.15 0.20 0.25
f(X) 1.0000 1.0513 1.1052 1.1618 1.2214 1.2840
∆2
∆ 513 539 566 596 626
26 27 30 30
∆3 1 3 0
Solución(a): f0′
=
1 0∆f0 h
−
1
∆2 f 0 2
= h 00.0513 −
+
∆3 f 0 1 3
0.00026 2
+
h = X i − X i−1 = 0.25 − 0.20 = 0.05
0.0001 1 3
f0′ = 1.0006667 ≅ 1.0007 Solución(b) ′ f(0.055) es una derivada que esta al principio de la tabla → se emplea la formula de Newton Progresivo, como XP=0.055 es un punto intermedio → XP ∈ [0.05,0.10+, → X0=0.05. Se halla X −X su valor de P de:P = P 0 . h
H=0.05 → P =
0.055−0.05 0.05
f0 X0 0.05 1.0513
Si fP′
=
f(′ 0.055 )
1 0∆f0 h
=
(2P−1)∆2 f 0 + 2!
= 0.1 ∆f0 539
∆2f0 27
∆3f0 3
(3P 2 −6P+2)∆2 f 0 + 1 3!
(0.2−1)(0.0027 ) 1 3(0.1)2 −6(0.1) +0 + 21 (0.0003)1 00.0539 + 0.05 2 6
= 1.05783
Ejemplo de la Fórmula Newton Regresivo Hallar(a) f ′ (0.025) y (b) f ′ (0.225) Solución (a) Para P=0 y ∇3 fi < 8, se utiliza → X 0
f0
0.25 1.289
′ fP=0
=
1 0∇f0 h
+
∇2 f 0 2
1
+
∇3 f 0 1 , f0′ 3
′ f(0.25) = f0′ = 0.05 00.0626 +
∇1 f0
∇2 f0
0.626
0.030
∇ 3 f0 0
′ ′ = fP=0 = f(0.25) , porque X0=0.25 en la tabla de NR
0.0030 2
+
0.0000 1 3
→ f0′ = 1.282
Solucion(b) X P −X 0 h
=
0.025−0.25 0.05
(2P+1)∇2 f 0 2!
+
(3P 2 +6P+2)∇3 f 0 1 3!
X P = 0.025 → P = 1
en fP′ = h 0∇f0 +
1
′ f(0.225) = 0.05 00.626 + 0 ′ f(0.225) = 1.252
= −0.5 → P = −0.5 → X0 f0 ∇1 f0 0.25 1.289 0.626
2(−0.5) 3(−0.5)2 +6(−0.5) (0.0030) + 0 + 11 + 21 (0.0000)1 2 2
∇ 2 f0 ∇ 3 f0 0.030 0
FÓRMULA DE DIFERENCIACIÓN DE BESSEL fP = f0 + β1 δf1/2 + β2 (δ2 f0 + δ2 f1 ) + β3 δ3 f1/2 + ⋯ 1 df P dP
fP′ = h
1
= h β1′ δf1/2 + β′2 (δ2 f0 + δ2 f1 ) + β′3 δ3 f1/2 + ⋯
Si β1 = P → β1′ P 2 −P 4
β2 =
P(P−1) 4
β3 =
P(P−1)(P−1/2) 6
=
→ β′2 =
=
2P−1 4
2P 3 −3P 2 +P 6
→ β′3 =
Si XP es un punto tabulado XP=X0 → P=0 1
f0′ = 0δf1/2 − h
(δ 2 f 0 +δ 2 f 1 ) 4
+
δ 3 f 1/2 4
+ ⋯1
6P 2 −6P+1 6
FÓRMULA DE DIFERENCIACIÓN DE STIRLING fp = f0 + s1 δf−1/2 + δf1/2 + s2 δ2 f0 + s3 δ3 f−1/2 + δ3 f1/2 + ⋯ fP′ =
1 h
s1′ δf−1/2 + δf1/2 + s2′ δ2 f0 + s3′ δ3 f−1/2 + δ3 f1/2 + ⋯ P
1
Si s1= 2 → s1′ = 2 , Si s2 =
P2 2
Si s3 =
→ s2′ = P ,
1 δf −1/2 +δf 1/2 1+ 2
fP′ = h 0
P(P−1)(P+2) 2(3!)
Si s4 =
Pδ2 f0 +
P 2 (P+1)(P−1) 4!
3P 2 −1 12
1 (δf −1/2 +δf 1/2 ) 2
−
(δ 3 f −1/2 +δ 3 f 1/2 ) 12
3P 2 12
→ s4′ =
2P 3 −P 12
δ3 f−1/2 + δ3 f1/2 + ⋯
En XP=X0 → P=0 f0′ = h 0
→ s3′ =
+ ⋯1
Ejemplo 3: Hallar(a) 𝑓 ′ (0.10), (b) 𝒇′ (𝟎. 125) Solución(a) f ′ (0.10) = f ′ (X0 ) = f0 Como 0.10 ∈ [0.10,0.15] → X 0 = 0.10 ∧ X p = 0.10 P=
X p −X 0 h
=
0.10−0.10 0.05
= 0 ∧ P = 0 ∈ < 0,1/2] se aplica f0′ de Stirling hasta su 3ra
diferencia porque segun T O T ∆3 f0 < 60 → f0′ de Stirling hasta su 3ra diferencia es: fP′ =
1 h
β1′ δf1/2 + β′2 (δ2 f0 + δ2 f1 ) + β′2 δ3 f1/2 … (i)
Como P = 0 → en (i)se tiene s3′ =
1
f0′ = h 0δf1/2 −
(δ 2 f 0 +δ 2 f 1 ) 4
+
δ 3 f 1/2 6
1
3P 2 −1 12
, s1′ =
1 2
,
s2′ = P
En la tabla se tiene: X 0 = 0.10
1.1052f0
X1 = 0.15
0.0𝟓𝟔𝟔δf1/2
0.00𝟐𝟕(δ2 f0 ) 0.0003(δ3 f1/2 ) 0.00𝟑𝟎(δ3 f1 )
1.1618f1
∧ h = 0.05 → f0′ =
1 0.0027+0.0030 − 00.0566 0.05 4
+
0.0003 1 6
= 1.1095
Solucion (b) f ′ (0.125) → X P = 0.125 ∈ [0.10,0.15] P=
X P −X 0 h
=
0.125−0.10 0.05
1
= 0.5 ∈ [2,1 > 𝑠𝑒 𝑎𝑝𝑙𝑖𝑐𝑎 𝐵𝑒𝑠𝑠𝑒𝑙 𝑎𝑠𝑡𝑎 𝑠𝑢 3𝑟𝑎 𝑑𝑖𝑓𝑒𝑟𝑒𝑛𝑐𝑖𝑎 𝑝𝑜𝑟
∆3 f0 < 60 𝑐𝑜𝑛 𝑃 = 0.5 𝑒𝑛 (i), donde ∶
β1′ = 1, β′2 =
2P −1 4
1
=
2(0.5)−1 4
= 0, β′3 =
6P 2 −6P+1 6
= −0.0833333 en (i)se tiene
f ′ (0.125) = 0.05 ,(1)(0.0566) + 0(0.0027 + 0.0030) + (−0.083)(0.0003)- = 1.1315
INTEGRACIÓN NUMERICA Para intervalos Simples Método del trapecio
=
𝑋0 − 𝑋1 𝑏−𝑎 = 1 1
𝑏
𝑓(𝑋)𝑑𝑥 = 𝑎
ℇ=−
3
(𝑓 + 𝑓1 ) + ℇ 2 0
𝑓 ′′ 12 (𝑚 )
,
m є *X0 , X1+
′′ ′′ ′′ Donde 𝑓(𝑚 ) = max { 𝑓(𝑎) , 𝑓(𝑏) }
Método de Simpson de 1/3
=
𝑏−𝑎 2
𝑏 𝑎
𝑓𝑥 𝑑𝑥 =
ℇ=−
5 90
(𝑓 + 4𝑓1 + 𝑓2 ) + ℇ 3 0
𝐼𝑉 𝑓(𝑚 )
,
m є [ X0 , X1 ]
𝐼𝑉 𝐼𝑉 𝐼𝑉 Donde 𝑓(𝑚 ) = max { 𝑓(𝑎) , 𝑓(𝑏) }
Método de Simpson de 3/8
= 𝑋3 𝑋0
𝑏−𝑎 3 3 𝑓(𝑥) 𝑑𝑥 = ,𝑓0 + 3𝑓1 + 3𝑓2 + 𝑓3 - + ℇ 8
ℇ=−
3 5 𝐼𝑉 𝑓(𝑚 ) 80
𝐼𝑉 𝐼𝑉 𝑑𝑜𝑛𝑑𝑒 𝑓(𝑚 ) = max { 𝑓(𝑎)
𝐼𝑉 , 𝑓(𝑏) }
EJEMPLO 2 𝑋 sin 𝑋 𝑑𝑥 1
=
2
( 𝑓0 + 𝑓1 ) + ℇ
Por el trapecio Simple
Solución:
X
𝑏−𝑎 2−1 = = =1 1 1
f(x)
X0=1
f0=0.841471
X0+h=X1=2
f1=1.818595
Obs: Los puntos Xi son dela forma Xi=X0+hi Determinación del ℇ Si 𝑓(𝑥) = 𝑋 sin 𝑋 ′ 𝑓(𝑥) = sin 𝑋 + 𝑋 cos 𝑋 ′′ 𝑓(𝑥) = 2 cos 𝑋 − 𝑋 sin 𝑋 ′′ ′′ 𝑓(′′𝑚 ) = max{ 𝑓(1) , 𝑓(2) }
𝑓(′′𝑚 ) = max{ 0.23 , 2.65 }
3
13
′′ ℇ = − 12 𝑓(𝑚 ) = − 12 (2.65) = −0.220907377
Entonces 2 𝑋 sin 𝑋 𝑑𝑥 1
=
2
( 𝑓0 + 𝑓1 ) + ℇ = 1.109125
Solución por el método Simpson de 1/3 para intervalo simple Sol:
=
𝑏−𝑎 2−1 1 = = = 0.5 2 2 2 X
f(x)
X0=1
f0=0.841471
X0+h=1.5
f1=1.496242
X0+2h=2
f2=1.818595
Determinación del ℇ 𝑓(′′𝑥 ) = 2 cos 𝑋 − 𝑋 sin 𝑋
Si
𝑓(𝐼𝑉 𝑥 ) = −4 cos 𝑋 + 𝑋 sin 𝑋 𝐼𝑉 𝐼𝑉 𝐼𝑉 𝑓(𝑚 { 𝑓(1) , 𝑓(2) } ) = max 𝐼𝑉 𝑓(𝑚 { 3.4831 } ) = max
ℇ=
5
𝐼𝑉 − 𝑓(𝑚 ) 90
=−
0.55 90
(3.4831) = −0.0012094
Entonces: 𝑏 𝑎
𝑓(𝑥) 𝑑𝑥 =
(𝑓 + 4𝑓1 + 𝑓2 ) + ℇ = 1.442048 3 0
Por Simpson 3/8 para intervalo Simple
Sol: 2 1
3 𝑋 sin 𝑋 𝑑𝑥 = (𝑓0 + 3𝑓1 + 3𝑓2 + 𝑓3 ) + ℇ 8
=
2−1 3
=
1 3
ℇ, es igual al ejemplo anterior
X
f(X)
X0=1
f0=0.841471
X1=4/3
f1=1.295917
X2=5/3
f2=1.659013
X3=2
f3=1.818595
𝐼𝑉 𝑓(𝑚 { 3.4831 } ) = max
ℇ=− 2 1
3 80
𝐼𝑉 5 𝑓(𝑚 ) = −
3 80
0.335 (3.4831) = −0.0005117
3 𝑋 sin 𝑋 𝑑𝑥 = (𝑓0 + 3𝑓1 + 3𝑓2 + 𝑓3 ) + ℇ = 1.42568923 8
Integración Numérica para intervalos compuestos
MÉTODO DEL TRAPECIO COMPUESTO
m veces h 𝑏−𝑎 = 𝑚 Demostración: 𝑋1 𝑋0 𝑋2 𝑋1 𝑋3 𝑋2
𝑓(𝑋)𝑑𝑥 =
(𝑓 + 𝑓1 ) + ℇ1 2 0
𝑓(𝑋)𝑑𝑥 =
(𝑓 + 𝑓2 ) + ℇ2 2 1
𝑓(𝑋)𝑑𝑥 =
(𝑓 + 𝑓3 ) + ℇ3 2 2 … …
𝑋𝑚
𝑓 (𝑋)𝑑𝑥 =
𝑋𝑚 −1
(𝑓 + 𝑓𝑚 ) + ℇ𝑚 2 𝑚 −1
Entonces la suma todas las integrales seria: 𝑋𝑚 𝑋0
𝑓(𝑋)𝑑𝑥 =
,(𝑓 + 𝑓𝑚 ) + 2(𝑓1 + 𝑓2 … … + 𝑓𝑚 −1 )- + ℇ 2 0
Determinación del ℇ del trapecio
(𝑏 − 𝑎)3 ′′ ℇ=− 𝑓 12𝑚2 (𝑚 )
Método de Simpson de 1/3 compuesta.
=
𝑏−𝑎 2𝑚
Demostración: 𝑋2 𝑋0 𝑋4 𝑋2 𝑋6 𝑋4
𝑋𝑚 −2
(𝑓 + 4𝑓1 + 𝑓2 ) + ℇ1 3 0
𝑓𝑥 𝑑𝑥 =
(𝑓 + 4𝑓1 + 𝑓2 ) + ℇ2 3 0
𝑓𝑥 𝑑𝑥 =
(𝑓 + 4𝑓1 + 𝑓2 ) + ℇ3 3 0
𝑓𝑥 𝑑𝑥 =
(𝑓 + 4𝑓1 + 𝑓2 ) + ℇ𝑚 3 0
…
𝑋𝑚
𝑓𝑥 𝑑𝑥 =
Entonces la suma de todas las integrales seria: 𝑿𝒎
𝒇(𝑿)𝒅𝒙 =
𝑿𝟎
𝒉 ,(𝒇𝟎 + 𝒇𝒎 ) + 𝟒(𝒇𝟏 + 𝒇𝟑 + 𝒇𝟓 + ⋯ + 𝒇𝒎−𝟏 ) + 𝟐(𝒇𝟐 + 𝒇𝟒 + ⋯ + 𝒇𝒎−𝟐 )- + ℇ 𝟑
Determinación del ℇ de Simpson de 1/3 compuesta
5 𝐼𝑉 ℇ1 = − 𝑓(𝑚 ) 90 5 𝐼𝑉 ℇ2 = − 𝑓(𝑚 ) 90
… … 5 𝐼𝑉 (𝑏 − 𝑎)5 𝐼𝑉 ℇ𝑚 = − 𝑓(𝑚 ) = 𝑓 90 2880𝑚4 (𝑚 )
Ejemplo: 2
Xsin 𝑋 𝑑𝑥 1
Por trapecio compuesto con m=5 Sol:
=
𝑏−𝑎 2−1 1 = = = 0.2 𝑚 5 5
ℇ=− 2 𝑓(𝑋)𝑑𝑥 1
=
2
(𝑏 −𝑎)3 12𝑚 2
′′ 𝑓(𝑚 ) =−
(2−1)3 12∗52
X
f(X)
X0=1
f0
X1=1.2
f1
X2=1.4
f2
X3=1.6
f3
X4=1.8
f4
X5=2
f5
′′ 𝑓(𝑚 )
,(𝑓0 + 𝑓5 ) + 2(𝑓1 + 𝑓2 + 𝑓3 + 𝑓4 )- + ℇ = 1.436070589 + ℇ
EXTRAPOLACION DE RICHARDSON- (E.R) Definición ER: Consiste en que a partir de dos estimación (o aproximaciones) de una integral, obtener una tercera aproximación (muchas veces la mejor). Se tiene los siguientes casos: Para intervalos simples: ER entre la Regla del trapecio y la 1º formula abierta (Integración abierta) “ambas de precisión uno” porque 𝑓 ′′ (ç) = 0 → 𝐸 = 0 Sea 𝐼 ∗ =
𝑏 𝑎
𝑓(𝑥)𝑑𝑥 por dos métodos anteriores se tiene:
Regla Trapecio (RT): I ∗ = I1 + E1 𝑑𝑜𝑛𝑑𝑒:
1
−3 = 12
( 1),
1
( , )
( 0,
1)
1
−( − )3 = 12
( 1),
1
1º Regla Abierta de Newton – Cotes: I ∗ = I2 + E2 𝑑𝑜𝑛𝑑𝑒:
2
−3 = 3
( 2),
2
( 0,
2)
1
= −
( , )
− 2
3
3
( 2),
El valor de la integral I* (generalmente es el mejor). Que se cumple: I ∗ = I1 + E1 = I2 + E2 … … … . (1) Tomando el cociente de errores (𝑏 − 𝑎)3 𝑖𝑖 𝑓 (ç1) 𝐸1 − 12 = (𝑏 − 𝑎)3 𝑖𝑖 𝐸2 (ç2) 24 𝑓
ç1 , ç2
(𝑎, 𝑏)
Supongamos que: f ii (ç1) = f ii (ç2) →
E1 24 = → E1 = −2E2 … … … . (2) E2 12
2
(2) en (1) I1 + E1 = I2 + E2 I1 − 2E2 = I2 + E2 I1 − I2 = E2 + 2E2 1
I1 − I2 = 3E2 → 3E2 = 3 (I1 − I2 ) (3) en (1) I ∗ = I 2 + E2 1
1
1
I ∗ = I2 + 3 (I1 − I2 ) = I2 + 3 I1 − 3 I2 1 3
1 3
I ∗ = I2 + I 2 ó
∗
=
1 3
1
+
2 … … … … … . (4) 3 2
E.R. entre la R. trapecio Simple y la primera formula abierta h =
b−a 2
EJEMPLO Calcular
3 3 (x 1
− 2x 2 + 7x − 5)dx = 1
20 3
→ F(x) = x 3 − 2x 2 + 7x − 5
2
Aplicando: 4 es decir: 𝐼∗ = 3 𝐼1 + 3 𝐼2 xi = x0 + 1h
i=1
i=2
x1 = x0 + 1h = 1 + 1(1) =2
x2 = x0 + 2h = 1 + 2(1) =3 h 2 I∗ = I1 + E1 = (f0 + f1 ) + E1 = (1 + 25) = 26 = I1 2 2
h=b−a =3−1= 2
X
𝑓(𝑥)
𝑋0 = 1
1 = 𝑓0
𝑋1 = 3
25 = 𝑓1
I = I2 + E2 = 2hf1 h=
b−a 2
=
3−1 2
f(X)
=1
f1
I = 2(1)(9) I2 = 18
Aplicando (4)
x0
1 2 I ∗ = I1 + I2 3 3 1
2
I ∗ = 3 (26) + 3 (18) =
62 3
x1
x2
2
= 20 3 X
𝑓(𝑥)
𝑋0 = 1
1 = 𝑓0
𝑋1 = 2
9 = 𝑓1
𝑋2 = 3
25 = 𝑓2
E.R. ENTRE LAS FORMULAS DE SIMPSON DE 1/3 Y 3/8 SIMPLE
Se sabe que para Simpson de 1/3 para intervalo simple es: I∗ = I1 + E1 =
x2 x0
5 𝐼𝑉 f(x)dx = (𝑓0 + 4𝑓1 + 𝑓2 ) − 𝑓 (£1), 3 90
b − a 5 1 IV → E1 = − f (£1), 2 90
£1
£1 (𝑥0 , 𝑥2 )
(b − a)5 IV (a, b) → E1 = − f (£1) 2880
Y que para Simpson de 3/8 es: ∗
I = I 2 + E2 =
x3 x0
3 35 𝐼𝑉 (𝑓 + 3𝑓1 + 3𝑓2 + 𝑓3 ) − f(x)dx = 𝑓 (£2), 8 0 80
3 b−a → E2 = − 80 3
5
(b − a)5 IV 1 IV (£2) f =− f (£2), £2 90 6480
(𝑎, 𝑏)
∴ 𝑰∗ = 𝑰𝟏 + 𝑬𝟏 = 𝑰𝟐 + 𝑬𝟐 … … (𝟏)` Sea el siguiente cociente de errores: (𝑏 − 𝑎)5 𝐼𝑉 (£1) − 𝐸1 2880 𝑓 = , £1 ≠ £2 𝑦 £1 ^ £2 (𝑏 − 𝑎)5 𝐼𝑉 𝐸2 − 6480 𝑓 (£2) Si 𝑓 𝐼𝑉 (£1) = 𝑓 𝐼𝑉 (£2) →
𝐸1 𝐸2
6480
(𝑎, 𝑏)
9
= 2880 → 𝐸1 = 4 𝐸2 … . . (2)`
9
9
(2)’ en (1)’: 𝐼1 + 4 𝐸2 = 𝐼2 + 𝐸2 → 𝐼2 − 𝐼1 = 4 𝐸2 − 𝐸2 →
5 4 𝐸2 = 𝐼2 − 𝐼1 → 𝐸2 = (𝐼2 −𝐼1 ) … … (3)` 4 5
£2
(𝑥0 , 𝑥3 )
(3)’ en (1)’: 4
I ∗ = I2 + E2 = I2 + 5 (I2 −I1 ) ∗
=
−
… (4)`
Formula De Extrapolación entre Simpson de 1/3 y Simpson de 3/8 para Intervalo Simple
Ejemplo: 9
4
Aplicando la fórmula I ∗ = 5 I2 − 5 I1 para hallar: 0.8
(0.2 + 25x − 200x 2 + 675x 3 − 900x 4 + 400x 5 )dx
0
Solución: 𝑃𝑎𝑟𝑎 𝐼1 : =
0.8 − 10 = 0.4 2
𝑋𝑖 = 𝑋0 + 𝑖 𝑋1 = 𝑋0 + 𝑖(0.4) F(x) = 0.2 + 25x − 200x 2 + 675x 3 − 900x 4 + 400x 5 I1 = I1 =
x2 f(x)dx x0
=
(𝑓 3 0
+ 4𝑓1 + 𝑓2 )
0.4 (0.2 + 4 ∗ (2.456) + 0.232) 3
I1 = 1.36746667
( ) 0
X0+1h =
1 2
0.4 0.8
0.2 2.456 0.232
0 1 2
( )
𝑋𝑖 = 𝑋0 + 𝑖
0 1
𝑋1 = 𝑋0 + 𝑖(0.4)
2 3
I2 = =
x3 f(x)dx x0 0.8−0 3
I2 =
=
x3 x0
0.8 3
=
3 8
0 0.2667 05333 0.8
0.2 1.43272438 3.48717696 0.232
(𝑓0 + 3𝑓1 + 3𝑓2 + 𝑓3 )
= 0.2667
f(x)dx =
3 0.2667(0.2 + 3(1.43272428) + 3(3.48717696) + 0.232) 8 9
𝐼2 = 1.51917037 →
9 5
4
𝐼 ∗ = 5 𝐼2 − 5 𝐼1 4 5
𝐼 ∗ = (1.51917037) − (1.36746667) 𝐼 ∗ = 164053333
0 1 2 3
(E.R) PARA EL TRAPECIO COMPUESTO
Sea:
𝐼 ∗ = 𝐼𝑛 1 + 𝐸𝑛 1 = 𝐼𝑛 2 + 𝐸𝑛 2 … … . (1)′′ Para 𝑛1 aplicaciones de la Regla Del Trapecio Compuesto, análogamente: 𝐼 ∗ = 𝐼𝑛 2 + 𝐸𝑛 2 Para 𝑛2 aplicaciones de la Regla Del Trapecio Compuesto tomando el cociente de errores: (𝑏 − 𝑎)3 𝐼𝑉 − 𝑓 (£1) 𝐸1 12𝑛22 = ; (𝑏 − 𝑎)3 𝐼𝑉 𝐸2 − 𝑓 (£2) 12𝑛12
𝑃𝑎𝑟𝑎 £1 ≠ £2 𝑦 £1 ^ £2
(𝑎, 𝑏)
Si:
𝑓
𝐼𝑉 (
£2) = 𝑓
𝐼𝑉 (
𝐸𝑛 2 𝑛1 £1) → = 𝐸𝑛 1 𝑛2
2
→ 𝐸𝑛 2
𝑛1 = 𝑛2
2
𝐸𝑛 1 … … . (2)′′
(2)” en (1)”: 2
𝑛
𝐼𝑛 1 + 𝐸𝑛 1 = 𝐼𝑛 2 + .𝑛 1 / 𝐸𝑛 1 2
𝑛 2 𝑛2
𝐼𝑛 2 − 𝐼𝑛 1 = 𝐸𝑛 1 − . 1 / 𝐸𝑛 1 𝑛
𝐸𝑛 1 1 − .𝑛 1 /
2
2
(3)” en (1)”:
𝐸𝑛 1 =
Para 𝑛2 = 2𝑛, en (3)”. Se tiene: 𝐸𝑛 1 =
𝐼𝑛 2 −𝐼𝑛 1
2
𝑛 1−. 1 /
→ 𝐸𝑛 1 =
2𝑛 1
= 𝐼𝑛 2 − 𝐼𝑛 1
𝐼𝑛 2 −𝐼𝑛 1 𝑛 1−. 1 /
2
… … . (3)′′
𝑛2
𝐼𝑛 2 −𝐼𝑛 1 3 4
4
𝐸𝑛 1 = 𝐼𝑛 2 − 𝐼𝑛 1 … … . (4)′′ 3
(4)” en (1)” 4
4
1
𝐼 ∗ = 𝐼𝑛 1 + 𝐸𝑛 1 = 𝐼𝑛 1 + 3 (𝐼𝑛 2 − 𝐼𝑛 1 ) = 3 𝐼𝑛 2 − 3 𝐼𝑛 1
∗
=
−
Extrapolación de Richardson (E.R.) para el Trapecio Compuesto para n1 y n2.
EJEMPLO Calcular
3 3 (𝑥 1
− 2𝑥 2 + 7𝑥 − 5)𝑑𝑥
usando “1º dos aplicaciones” (𝑛1 = 2) y despues “cuatro
aplicaciones” (𝑛2 = 4) de la regla del Trapecio compuesto y despues aplicar E.R. para encontrar una 3º aproximacion (que es la mejor).
Solución:
𝐼𝑛 = ,𝑓0 + 𝑓𝑛 + 2(𝑓1 + 𝑓2 + 𝑓3 + ⋯ + 𝑓𝑛−1 - + (𝐸 ) 2
Se tiene: 𝑛1 = 2 y 𝑛2 = 4 Para:
𝐼𝑛 = 2 (𝑓0 + 𝑓𝑛 + 2
𝑛−1 𝑖=1 𝑓𝑖 )
Como 𝑛1 = 2
→ =
𝑏−𝑎 𝑛1
Para:
1 9 2.5
3−1 2
=1 1
𝐼𝑛 1 = 2 (1 + 25 + 2(9))
( ) 1 2 3
=
0 1 2
𝐼𝑛 1 =
44 2
= 22
𝐼𝑛 = 𝑓 + 𝑓𝑛 + 2 2 0
𝑛−1
𝑓𝑖 𝑖=1
Como 𝑛2 = 4
→ =
𝑏−𝑎 𝑛2
=
3−1 4
= 0.5
( ) 1 1.5 2 2.5 3
1 4 9 15 25
xi = x0 + ih
0 1
x1 = 1 + 1(0.5)
2 3 4
x1 = x0 + 1h
𝐼𝑛 2 =
0.5 2
3
35
𝐼𝑛 2 = 0.25 .26 + 2( 8 + 𝐼𝑛 2 = 21
5
.1 + 25 + 2(4 8 + 9 + 15 8)/ 125 8
)/
Aplicando: 𝟒
𝟏
𝐈 ∗ = 𝟑 𝐈𝐧𝟐 − 𝟑 𝐈𝐧𝟏 𝟒
𝟏
𝐈 ∗ = 𝟑 (𝟐𝟏) − 𝟑 (𝟐𝟐) = 𝐈 ∗ = 𝟐𝟎
𝟐 𝟑
𝟔𝟐 𝟑
INTEGRACION DE ROMBERG Es un método que combina la convergencia de las reglas del Trapecio o Simpson con las técnicas de extrapolación de Richardson para retener una sucesión que converge hacia el verdadero valor de la integral. Descripción de la técnica de integración de Romberg Vamos a introducir la 𝑅𝐾 , la aproximacion de la integral 𝐼 = Trapecio Compuesto para 𝑛𝐾 = 2𝑘−1 y 𝐾 = Para: 𝑛𝐾 = 2𝑘−1 y 𝐾 =
𝑏−𝑎 𝑛𝑘
=
𝑏−𝑎 2𝑘 −2
𝑏−𝑎 𝑛𝑘
=
𝑏−𝑎 2𝑘 −2
𝑏 𝑎
𝑓(𝑥)𝑑𝑥 utilizando la regla del
.
; se sugiere la siguiente tabla:
𝑲
1
2
3
4
𝒏𝑲 = 𝟐𝑲−𝟏
1
2
4
8
𝑏−𝑎 1
𝑏−𝑎 2
𝑏−𝑎 4
𝑏−𝑎 8
𝒉𝑲 =
𝒃−𝒂 𝒏𝒌
La integración de Romberg da su respuesta en forma matricial:
𝑅11 𝑅21
𝑅22
𝑅31
𝑅32
𝑅33
𝑅41
𝑅42
𝑅43
𝑅44
𝑅51
𝑅52
𝑅53
𝑅54
. . .
𝑅𝑛1
. . .
𝑅𝑛2
. . .
𝑅𝑛3
. . .
𝑅𝑛4
𝑅55 . . .
𝑅𝑛5
.
. .
𝑅𝑛𝑛
…
n
Se sugiere hallar 𝑅11 con la siguiente formula.
1
= 2
( )+
2(
( )+ 2
)
1
( + )
Nos generamos las siguientes sucesiones 𝑅11 , 𝑅21 , 𝑅31 , … 𝑅𝐾1
𝑘=1
𝑥1 𝑥0
𝑓 (𝑥 )𝑑𝑥 = 𝑅11 =
1 (𝑓(𝑎) − 𝑓(𝑏)) 2
Obs: Para las demás filas y columnas se sugiere hallar con la siguiente formula
1
1 = 2
( −1),1
+
2 −1
1
1 ( + ( − ) 2
−1 )
De la segunda columna hacia adelante se calculara con la siguiente fórmula:
=
4
−1
, −1 − 4 −1 −
−1, −1
1
i= 2,3,4,…n j=
Para la segunda columna se utilizara: 2
=
4
− 3
1
−1,1
k = 2,3,…
Para la tercera columna k=3 será: 3
=
16
2
− 15
−1,2
k = 3,4,…
Para la cuarta columna k=4 será: 4
=
64
3
− 63
−1,3
k = 4,5,…
Para la quinta columna k=5 será: 5
=
256
− 255 4
−1,4
k = 5,6,…
Para la sexta columna k=6 será: 6
=
1024
− 1023 5
−1,5
k = 6,7,…
EJEMPLO Hallar
𝟎.𝟖 𝒙𝒆𝒙 𝒅𝒙 ; 𝟎
con n=6 por Romberg
Solución: Tenemos de datos: a=0 b=0.8 n=6 𝑓(𝑥) = 𝑥𝑒 𝑥 Realizando la tabla 𝑲
1
2
3
4
5
6
𝒏𝑲 = 𝟐𝑲−𝟏
1
2
4
8
16
32
0.8
0.4
0.2
0.1
0.05
0.025
𝒉𝑲 =
𝒃−𝒂 𝒏𝒌
Hallamos la primera fila y columna
𝑅1𝐾
𝐾 = 𝑓 (𝑎 ) + 𝑓 (𝑏 ) + 2 2
2𝑘−1 𝑖=1
𝑓(𝑎 + 𝑖𝐾 )
Cuando K = 1
𝑅11 =
1
𝑅11 =
1
2
2
𝑓 (0) + 𝑓(0.8) + 2
21−1 𝑖=1
𝑓(𝑎 + 𝑖1 )
,𝑓 (0) + 𝑓 (0.8)-
𝑓 (𝑥) = 𝑥𝑒 𝑥 𝑓(0) = 0𝑒 0 = 0 𝑓(0.8) = 0.8𝑒 0.8 = 1.780432743
𝑅11 =
0.8 2
,0 + 1.780432743- = 0.7121730971
𝑅11 = 0.7121730971
Calculando la primera columna y demás filas. 𝑅𝐾1
1 = 𝑅( ) + 𝐾−1 2 𝐾−1 ,1
2𝐾−2 𝑖=1
1 𝑓(𝑎 + (𝑖 − )𝐾−1 ) 2
Cuando K = 2
𝑅21 =
1
0𝑅(2−1),1 + 2−1 2
22−2 𝑖=1
𝑓(𝑎 + (𝑖 − 2)2−1 )1 = 0.594777
1
1
0𝑅(3−1),1 + 3−1 2
23−2 𝑖=1
𝑓(𝑎 + (𝑖 − 2)3−1 )1 = 0.564899
1
0𝑅(4−1),1 + 4−1 2
24−2 𝑖=1
𝑓(𝑎 + (𝑖 − 2)4−1 )1 = 0.557395
1
0𝑅(5−1),1 + 5−1
25−2 𝑖=1
𝑓(𝑎 + (𝑖 − )5−1 )1 = 0.555517
0𝑅(6−1),1 + 6−1
26−2 𝑖=1
𝑓(𝑎 + (𝑖 − )6−1 )1 = 0.555047
Cuando K = 3
𝑅31 =
1
Cuando K = 4
𝑅41 =
1
Cuando K = 5
𝑅51 =
2
1 2
Cuando K = 6
𝑅61 =
1 2
1 2
PARA LOS VALORES DE LA SEGUNDA COLUMNA:
𝑅𝑘2 =
4𝑅𝑘1 −𝑅𝑘−1,1 3
Cuando K = 2
𝑅22 =
4𝑅21 −𝑅2−1,1 3
= 0.555545
Cuando K = 3
𝑅32 =
4𝑅31 −𝑅3−1,1 3
= 0.554939
Cuando K = 4
𝑅42 =
4𝑅41 −𝑅4−1,1 3
= 0.554893
Cuando K = 5
𝑅52 =
4𝑅51 −𝑅5−1,1 3
= 0.554891
Cuando K = 6
𝑅62 =
4𝑅61 −𝑅6−1,1 3
= 0.554890
PARA LOS VALORES DE LA TERCERA COLUMNA:
𝑅𝑘3 =
16𝑅𝑘2 −𝑅𝑘−1,2 15
Cuando K = 3
𝑅33 =
16𝑅32 −𝑅3−1,2 15
= 0.554891
Cuando K = 4 𝑅43 =
16𝑅42 −𝑅4−1,2 15
= 0.554889
Cuando K = 5
𝑅53 =
16𝑅52 −𝑅5−1,2 15
= 0.554890
Cuando K = 6
𝑅63 =
16𝑅62 −𝑅6−1,2 15
= 0.554889
PARA LOS VALORES DE LA CUARTA COLUMNA:
𝑅𝑘4 =
64𝑅𝑘3 −𝑅𝑘−1,3 63
Cuando K = 4
𝑅44 =
64𝑅43 −𝑅4−1,3 63
= 0.554888
Cuando K = 5
𝑅54 =
64𝑅53 −𝑅5−1,3 63
= 0.554890
Cuando K = 6
𝑅64 =
64𝑅63 −𝑅6−1,3 63
= 0.554889
PARA LOS VALORES DE LA QUINTA COLUMNA:
𝑅𝑘5 =
256𝑅𝑘4 −𝑅𝑘−1,4 255
Cuando K = 5
𝑅55 =
256𝑅54 −𝑅5−1,4 255
= 0.554890
Cuando K = 6
𝑅65 =
256𝑅64 −𝑅6−1,4 255
= 0.554889
PARA LOS VALORES DE LA SEXTA COLUMNA:
𝑅𝑘6 =
1024 𝑅𝑘5 −𝑅𝑘−1,5 1023
Cuando K = 6
𝑅66 =
1024 𝑅65 −𝑅6−1,5 1023
= 0.554889
LA MATRIZ FINAL SERÁ:
0.712173 0.594777
0.55545
0.564899
0.554939
0.554891
0.557395
0.554889
0.554889
0.554888
0.555517
0.554890
0.554890
0.554890
0.554890
0.555047
0.554890
0.554889
0.554889
0.554889
0.554889
SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES ORDINARIAS CON VALOR INICIAL.
I) De un solo paso: Se tiene los siguientes métodos: a) Método de Euler b) Método de Taylor c) Método de Runge-Kutta de orden 4 d) Método de Runge-Kutta de orden 2 *Todas ellas son de la siguiente forma: Dada la Ecuación Diferencial Ordinaria:
= ( , ) /
=
Con la condición inicial:
( 0) =
0
a) Método de Euler: Es de la forma: +1
=
+
O
+1
=
+ ( , )
Observación: Euler es un caso particular de Taylor de orden 1.
ejemplo 1º forma de pregunta 1) Dado: 𝑦´ = 2𝑥𝑦 2 Con: 𝑦(1) = 1 Hallar: 𝑦(1.1)
Solución: Como: 𝑥0 = 1
, 𝑦0 = 1
Entonces: Para i=0 , 𝑦1 = 𝑦0 + 𝑦´0 Donde: 𝑦1 = 𝑦(1.1) = 𝑦(𝑥1 ) → 𝑥1 = 1.1 → = 𝑥1 − 𝑥0 = 1.1 − 1 = 0.1 = 0.1 Luego: 𝑦1 = 𝑦0 + 𝑦´0 … (∝)
Donde: 𝑦´0 = 2𝑥0 𝑦02 Reemplazando en (∝): 𝑦1 = 𝑦0 + 2𝑥0 𝑦02 𝑦1 = 1 + (0.1)(2)(1)(1)2 𝑦1 = 1.2 𝑦(1.1) = 1.2
2º forma de pregunta 2) Dado: 𝑦´ = 2𝑥𝑦 2 Con: 𝑦(1) = 1 Hallar: 𝑦(1.1) para 𝑥 ∈ [1, 2] , con n=4 Solución: =
Como:
𝑏−𝑎 𝑛
=
2−1 4
= 0.25
𝑖
0
1
2
𝑥𝑖
𝑥0
𝑥1
𝑥2
𝑥3
1.25 1.5
1.75
1
3
4 𝑥4 2
𝑥1 = 𝑥0 + 1(0.25) = 1.25 𝑥2 = 𝑥0 + 2(0.25) = 1.5 𝑥3 = 𝑥0 + 3(0.25) = 1.75 𝑥4 = 𝑥0 + 4(0.25) = 2
Para i=0 , con (𝑥0 , 𝑦0 ) = (1, 1) 𝑦1 = 𝑦0 + 𝑦´0 𝑦1 = 𝑦0 + 2𝑥0 𝑦02 𝑦1 = 1.2
Para i=1 , con (𝑥1 , 𝑦1 ) = (1.25, 1.2) 𝑦2 = 𝑦1 + 2𝑥1 𝑦12 𝑦2 = 2.1
Para i=2 , con (𝑥2 , 𝑦2 ) = (1.5, 2.1) 𝑦3 = 𝑦2 + 2𝑥2 𝑦22 𝑦3 = 5.4075 Para i=3 , con (𝑥3 , 𝑦3 ) = (1.75, 5.4075) 𝑦4 = 𝑦3 + 2𝑥3 𝑦32 𝑦4 = 30.99342422
Para i=4 , con (𝑥4 , 𝑦4 ) = (2, 30.99342422) 𝑦5 = 𝑦4 + 2𝑥4 𝑦42 𝑦5 = 991.5857691
b) Método de Taylor de orden K:
Es de la siguiente forma: (𝑘)
𝑦𝑖+1
𝑘 𝑦𝑖 𝑦´𝑖 2 𝑦𝑖´´ 3 𝑦𝑖´´´ = 𝑦𝑖 + + + +⋯ + 1! 2! 3! 𝑘!
+ 𝜀
O 𝑦 ´ (𝑥𝑖 ) 2 𝑦 ´´ (𝑥𝑖 ) 3 𝑦 ´´´ (𝑥𝑖 ) 𝑘 𝑦 (𝑘) (𝑥𝑖 ) 𝑦(𝑥𝑖+1 ) = 𝑦(𝑥𝑖 ) + + + +⋯ + + 𝜀 1! 2! 3! 𝑘!
Donde: El error local es de la forma: 𝑘+1 𝑦 (𝑘+1) (𝑥𝑖 ) 𝜀= 𝑘 + 1!
1º forma de pregunta Dado: 𝑦 ´ = 𝑥 2 + 𝑦 2 … (∝) Con: 𝑦(1) = 0 Por Taylor de orden 3 hallar: 𝑦(1.5)
Solución: Como: 𝑥0 = 1 𝑦 𝑦0 = 0 Y de 𝑦(1.5) se sabe que: 𝑥1 = 1.5 → = 𝑥1 − 𝑥0 = 1.5 − 1 = 0.5 Luego, Taylor de orden 3 es de la forma: 𝑦𝑖+1 Para i=0 con (𝑥0 , 𝑦0 ) = (1, 0) 𝑦1 = 𝑦0 +
𝑦´0 1!
+
2 𝑦0´´ 2!
Se requiere hallar:
+
3 𝑦0´´´ 3!
𝑦´𝑖 2 𝑦𝑖´´ 3 𝑦𝑖´´´ = 𝑦𝑖 + + + 1! 2! 3!
𝑦´0 , 𝑦0´´ , 𝑦0´´´ 𝐷𝑒 (∝): 𝑦 ´ = 𝑥 2 + 𝑦 2 … (𝛽) 𝑦0´ = 𝑥02 + 𝑦02 𝑦0´ = 12 + 02 𝑦0´ = 1 Derivando (𝛽): 𝑦 ´´ = 2𝑥 + 2𝑦𝑦 ´ … (𝜃) 𝑦0´´ = 2𝑥0 + 2𝑦0 𝑦0´ 𝑦0´´ = 2(1) + 2(0)(1) 𝑦0´´ = 2 Derivando (𝜃): 𝑦 ´´´ = 2 + 2(𝑦𝑦 ´´ + 𝑦 ´ 𝑦 ´ ) 𝑦0´´´ = 2 + 2(𝑦0 𝑦0´´ + 𝑦0´ 𝑦0´ ) 𝑦0´´´ = 2 + 2((0)(2) + (1)(1)) 𝑦0´´´ = 4
Por último, reemplazando: 𝑦1 = 0 +
80.5)(0) (0.5)2 (2) + 1! 2!
+
(0.5)3 (4) 3!
𝑦1 = 0.8333333 … 𝑦(1.5) = 0.8333333 …
2º forma de pregunta Dado: 𝑦 ´ = 𝑥 2 + 𝑦 2 … (∝) Con: 𝑦(1) = 0 Por Taylor de orden 3 en el segmento 𝑥 ∈ [1, 2], con n=2, hallar: 𝑦(1.5) Solución: Como:
=
𝑏−𝑎 𝑛
=
2−1 2
= 0.5
Determinamos los puntos: 𝑥𝑖 = 𝑥0 + 𝑖 / 𝑥𝑖 ∈ ,1, 2-
𝑖 = 0 → 𝑥0 = 𝑥0 + 0 = 1 𝑖 = 1 → 𝑥1 = 𝑥0 + 1 = 1.5 𝑖 = 2 → 𝑥2 = 𝑥0 + 2 = 2
Para i=0 con (𝑥0 , 𝑦0 ) = (1, 0) 𝑦1 = 𝑦0 +
𝑦´0 1!
+
2 𝑦0´´ 2!
𝑦1 = 0.8333333 …
+
3 𝑦0´´´ 3!
(Resuelto anteriormente)
Para i=1 con (𝑥1 , 𝑦1 ) = (1.5, 0.8333333)
𝑦2 = 𝑦1 +
𝑦´1 1!
+
´´ 2 𝑦1 2!
+
´´´ 3 𝑦1 3!
Se requiere hallar: 𝑦´1 , 𝑦1´´ , 𝑦1´´´ De (∝): 𝑦 ´ = 𝑥 2 + 𝑦 2 …(𝛽) 𝑦1´ = 𝑥12 + 𝑦12 𝑦1´ = 1.52 + 0.8333333 …2 𝑦1´ = 2.9444444 … Derivando (𝛽): 𝑦 ´´ = 2𝑥 + 2𝑦𝑦 ´ … (𝜃) 𝑦1´´ = 2𝑥1 + 2𝑦1 𝑦1´ 𝑦1´´ = 2(1.5) + 2(0.8333333 … )(2.9444444 … ) 𝑦1´´ = 7.907407405 Derivando (𝜃): 𝑦 ´´´ = 2 + 2(𝑦𝑦 ´´ + 𝑦 ´ 𝑦 ´ ) 𝑦1´´´ = 2 + 2(𝑦1 𝑦1´´ + 𝑦1´ 𝑦1´ ) 𝑦1´´´ = 2 + 2((0.8333333 … )(7.907407405) + (2.9444444 … )(2.9444444 … )) 𝑦1´´´ = 32.5185185
Luego , reemplazando: 𝑦2 = 𝑦1 +
𝑦´1 1!
+
2 𝑦1´´ 2!
𝑦2 = 0.8333333 … +
+
3 𝑦1´´´ 3!
(0.5)(2.9444444 … ) 1!
+
(0.5)2 (7.907407405 )
𝑦2 = 3.235339504
Para i=2 con (𝑥2 , 𝑦2 ) = (2, 3.235339504 ) 𝑦3 = 𝑦2 +
𝑦´2 1!
+
2 𝑦2´´ 2!
+
3 𝑦2´´´ 3!
Se requiere hallar: 𝑦´2 , 𝑦2´´ , 𝑦2´´´ De (∝):
2!
+
(0.5)3 (32.5185185 ) 3!
𝑦 ´ = 𝑥 2 + 𝑦 2 …(𝛽) 𝑦2´ = 𝑥22 + 𝑦22 𝑦2´ = 22 + 3.2353395042 𝑦2´ = 14.46742171 Derivando (𝛽): 𝑦 ´´ = 2𝑥 + 2𝑦𝑦 ´ … (𝜃) 𝑦2´´ = 2𝑥2 + 2𝑦2 𝑦2´ 𝑦2´´ = 2(2) + 2(3.235339504)(14.46742171) 𝑦2´´ = 97.61404193
Derivando (𝜃): 𝑦 ´´´ = 2 + 2(𝑦𝑦 ´´ + 𝑦 ´ 𝑦 ´ ) 𝑦2´´´ = 2 + 2(𝑦2 𝑦2´´ + 𝑦2´ 𝑦2´ ) 𝑦2´´´ = 2 + 2((3.235339504)(97.61404193) + (14.46742171)(14.46742171)) 𝑦2´´´ = 1052.241714
Luego , reemplazando: 𝑦3 = 𝑦2 +
𝑦´2 1!
+
2 𝑦2´´ 2!
𝑦3 = 3.235339504 + 𝑦3 = 44.59250797
+
3 𝑦2´´´ 3!
(0.5)(14.46742171 ) 1!
+
(0.5)2 (97.61404193 ) 2!
+
(0.5)3 (1052 .241714 ) 3!
PREDICTOR – CORRECTOR Como su nombre lo indica 1𝑟𝑜 se predice un valor 𝑦𝑚+1 , después se usa una formula diferente para corregir este valor. Los métodos Predictor-Corrector hallar el valor del pto (𝑥𝑚+1 , 𝑦𝑚+1 ) utiliza información del pto (𝑥𝑚 , 𝑦𝑚 ) y sus precedentes donde cada pto precedente indica un paso, de aquí su nombre como método de múltiples pasos . A diferencia de los métodos de Runge-Kuta, Taylor, Euler, que son métodos de un solo paso (Para 𝑦𝑖+1 utiliza la información de 𝑦𝑖 )
ejemplo Usando el método predictor – corrector
Predictor: Corrector:
= +1 =
+1
+ 2 ′ + ( ′+ 2
−1
′ +1 )
Con h=0.5, hallar y(1) sabiendo que y satisface: y ′ = xy 2 + 3x , 𝑦0 = −0.5 Solución: y ′ = f(x,y) = xy 2 + 3xy ;
𝑦0 = −0.5 Y h=0.5
Por condición del problema. 𝑦𝑖−1 = 𝑦0 = −0.5 𝑦𝑖 = 𝑦(0.5) 𝑦𝑖+1 = 𝑦(1) Para aplicar el predictor corrector, se necesita hallar 𝑦1 por un método no menor de 2𝑑𝑜 orden (puede ser Taylor, o RK), Sea RK-4.
𝑦1 = 𝑦0 + (𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 ) 𝑏 ′ 𝑘1 = 𝑓(𝑥 0 ,𝑦0 ) = y0′ = y(0) =0 𝑒𝑛 y ′ = xy 2 + 3xy y0′ = 𝑥0 y02 + 3𝑥0 𝑦0 = 0(−0.5)2 + 3(0)(−0.5) y0′ = 0 𝑘2=𝑓(𝑥
0 +2
1 2
, 𝑦0 + 𝑘 1 )
= 𝑓(0.25,00.5+0.25(0))
𝑘2 = 𝑓(0.25,−0.5) = (0.25)(−0.5)2 + 3(0.25)(−0.5) = −0.3125 𝑘3 = 𝑓
1 (𝑥 0 + , 𝑦0 + 𝑘 2 ) 2 2
= 𝑓(0.25 ,−0.5+(0.25)(−0.3125))
𝑘3 = 𝑓(0.25,−0.578125 ) 𝑘3 = 0.25(−0.578125)2 + 3(0.25)(−0.578125) = −0.3500366
𝑘4 = 𝑓
1 (𝑥 0 + , 𝑦0 + 𝑘 3 ) 2 2
𝑘4 = 𝑓(0.5,−0.5+0.5(−0.350066 )) 𝑘4 = 𝑓(0.5,−0.67533 ) 𝑘4 = 0.5(−0.67533)2 + 3(0.5)(−0.67533) = −0.7847026
𝑆𝑖 𝑦1 = 𝑦0 + 6 (𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 )
0.5 𝑦1 = −0.5 + (… ) 6 𝑦𝑖 = 𝑦1 = −0.6758146 𝑦𝑖′ = (0.5)(−0.6758146)2 − 3(0.5)(−0.6758146) 𝑦𝑖′ = −0.7853592
Aplicando el predictor: 𝑦𝑖+1 = 𝑦𝑖−1 + 2𝑦𝑖′ 𝑦(1) = −0.5 + 2(0.5)(−0.7853592) 𝑦(1) = −1.2853592 ′ ′ 𝑦𝑖+1 = 𝑦(1) = 1(−1.2853592)2 + 3(1)(−1.2853592) ′ 𝑦𝑖+1 = −2.2039293
Aplicando el corrector: ′ 𝑦𝑖+1 = 𝑦𝑖 + (𝑦𝑖′ + 𝑦𝑖+1 ) 2 0.5 𝑦(1) = −0.6758146 + (−0.7853592 ± 2.2039293) 2 𝑦(1) = −1.423167 → este es el corrector
ECUACIONES DIFERENCIALES DE ORDEN SUPERIOR Una ecuación diferencial ordinaria de orden n se puede transformar en un sistema de EDO de orden 1. Así se tiene una EDO de orden n:
( )
=
( , ,
(𝛿)
,
,
,…, (
)
)
𝛿 con las condiciones iniciales 𝑦(𝑥 = 𝑦0 n condiciones iniciales 0)
Sea las sgtes transformaciones 𝑦1 = y 𝑦2 = y ′ 𝑦3 = y ′′ . 𝑦n = y (n−1) Se tiene 𝑦1′ = y ′ = 𝑦2 𝑦2′ = y ′′ = 𝑦3 . 𝑦𝑛′ = 𝑦n = 𝑓(𝑥,𝑦 ,𝑦 ′ ,𝑦 ′′ ,𝑦 ′′′ ,…,𝑦 (𝑛 −1) )
Sea 𝑦i = 𝑦(𝑥 i ) =
𝑦1(𝑥 i ) 𝑦2(𝑥 i ) … 𝑦n(𝑥 i )
′ → 𝑦𝑖′ = 𝑦(𝑥 = 𝑖)
′ 𝑦1(𝑥 i) ′ 𝑦2(𝑥 i ) … ′ 𝑦𝑛(𝑥 i)
Ejemplo 1: La Sgte. EDO de 3𝑒𝑟 orden 𝑦 ′′′ + 𝑡𝑦 ′′ − 𝑡𝑦 ′ = 𝑡 ′′ ′ Con 𝑦(0) = 𝑦(0) = 0 , 𝑦(0) =1
Como un conjunto de ecuaciones de 1 𝑒𝑟 orden Solución: De 𝑦 ′′′ + 𝑡𝑦 ′′ − 𝑡𝑦 ′ − 2𝑦 = 𝑡 → 𝑦 ′′′ = 𝑡 + 2𝑦 + 𝑡y ′ − 𝑡y ′′ Sea 𝑦 = 𝑦1 𝑦1(0) = 𝑦(0) = 0 ′ 𝑦1′ = 𝑦 ′ = 𝑦2 𝑦2(0) = 𝑦(0) =1 𝑦2′ = 𝑦 ′′ = 𝑦3 𝑦3′ = 𝑦 ′′′
′′ 𝑦3(0) = 𝑦(0) =0
La EDO de3𝑒𝑟 orden se expresa como un conjunto de 3 ecuaciones de 1𝑒𝑟 orden sgte: 𝑦1′ = 𝑦2 𝑦2′ = 𝑦3 𝑦3′ = 𝑡 + 2𝑦1 + 𝑡𝑦2 − 𝑡𝑦3
Matricialmente se tiene: = 𝑦′ = 𝑦1′ 𝑦2 ′ 𝑦2 𝑦3 ′ 𝑦3 𝑡 + 2𝑦1 + 𝑡𝑦2 − 𝑡𝑦3 Con las Sgts. condiciones iniciales 𝑦(0) =
𝑦1(0) 𝑦2(0) 𝑦3(0)
=
0 1 0
