• Author / Uploaded
• Oscar Angulo

Citation preview

MÉTODOS NUMÉRICOS Y DISEÑO DE ALGORITMOS I

UNIVERSIDAD MAYOR DE SAN ANDRÉS FACULTAD DE INGENIERÍA INGENIERÍA MECÁNICA-ELECTROMECÁNICA

Métodos Numéricos INGENIERÍA MECÁNICA-ELECTROMECÁNICA FACULTAD DE INGENIERÍA y UNIVERSIDAD MAYOR DE SAN ANDRÉS Diseño de Algoritmos PROBLEMAS RESUELTOS

MEC-130 Docente: ING. PASTOR L. BARRÓN L. Autor: UNIV. MARTIN LAGUNA A [email protected]

La Paz, Bolivia

Julio 2019

MÉTODOS NUMÉRICOS MEC-130

Martin Laguna A. Universidad Mayor de San Andrés Facultad de Ingeniería Ingeniería Mecánica-Electromecánica Julio 2019

Tabla de Contenido

2

Capítulo 1 Primer Parcial 1.1 1.2 1.3

17

Segundo Parcial SERIES DE FOURIER 17 MÍNIMOS CUADRADOS 21 AJUSTE MULTILINEAL 23

Capítulo 3 Recuperatorio 3.1 3.2

33

10

Capítulo 2

2.1 2.2 2.3

26

NEWTON RAPHSON 2 GAUSS SEIDEL 6 ELEMENTOS FINITOS EN UNA DIMENSIÓN

POLINOMIO DE LAGRANGE 26 POLINOMIO DE NEWTON (DIFERENCIAS DIVIDIDAS)

Capítulo 4 Examen Final 4.1 4.2 4.3

ROMBERG 33 CUADRATURA DE GAUSS CUADRATURA DE GAUSS

39 43

30

2

1 1.1

Primer Parcial

NEWTON RAPHSON

Ejercicio 1 La ecuación de una onda estacionaria reflejada esta dada por: ¶ µ ¶ ¸ · µ 2πt v 2πx h = h 0 sen cos + e −x λ λ

Donde h = 0,4 · h0 , λ = 16, t = 12, v = 48 , calcular el valor más pequeño de x .

Desarrollo y solución

E

l problema se reduce a encontrar las raíces de la ecuación de la forma f (x) = 0.

La ecuación también puede escribirse como: · µ ¶ µ ¶ ¸ h 2πx 2πt v = sen cos + e −x h0 λ λ ·

µ sen

¶ µ ¶ ¸ 2πx 2πt v h cos + e −x − =0 λ λ h0

(1.1)

entonces, para una h, h0 , λ, t y v la función en variable x viene dada por la ecuación: ¶ µ ¶ ¸ 2πx 2πt v h −x f (x) = sen cos +e − λ λ h0 ·

µ

(1.2)

reemplazando los valor numéricos se obtiene. f (x) = sen (0,392699x) + e −x − 0,4

Auxiliar: Martin Laguna A.

(1.3)

3 1.1. NEWTON RAPHSON

N

para encontrar la raíz vamos a evaluar la ecuación (1.3) en un intervalo de confianza.

inicio = 0

incremento = 1

;

;

fin = 10

para así obtener la siguiente tabla. i

xi

y i = f (x i )

1

0

0.6

2

1

0.35056

3

2

0.44244

4

3

0.57367

5

4

0.61832

6

5

0.53062

7

6

0.30959

8

7

-0.016405

9

8

-0.39966

10

9

-0.78256

11

10

-1.1071

Tabla. 1.1

de la Tabla 1.1 se puede observar que existe un cambio de signo en el intervalo [x = 6, x = 7], por tanto una raíz se encuentra entre estos valores. Para resolver el problema tomaremos los siguientes parámetros.

x0

Decimales

Tolerancia

7

5

10−3

Tabla. 1.2

Planteamos la relación de recurrencia. f (x n ) ˙ n) f (x sen (0,392699x n ) + e −xn − 0,4 x n+1 = x n − 0,392699cos (0,392699x n ) − e −xn x n+1 = x n −

Por tanto con la tabla 1.2 y la ecuación (1.4) procedemos a realizar las iteraciones.

Primera Iteración:

Auxiliar: Martin Laguna A.

(1.4)

4 1.1. NEWTON RAPHSON

Para n = 0 determinamos x 1 .

x1

= x0 −

sen (0,392699x 0 ) + e −x0 − 0,4 0,392699cos (0,392699x 0 ) − e −x0

= 7−

sen (0,392699 · 7) + e −7 − 0,4 0,392699cos (0,392699 · 7) − e −7

=

6,954897

por tanto el error

Er r or 1

= |6,954897 − 7| =

0,045103

Error1 |−5| ; |21| > |−2|

|22| > |−5| ; |22| > |−5|

A partir del sistema de ecuaciones lineales, se debe despejar a la incógnita que se ubica en la diagonal principal de cada una de las ecuaciones que conforma el sistema y así obtenemos las relaciones de recurrencia que nos permitan realizar las iteraciones correspondientes para determinar la solución del mismo: x (k+1) =

2y (k) 3z (k) 500 + + 17 17 17

y (k+1) =

5x (k) 2z (k) 200 + + 21 21 21

z (k+1) =

5x (k) 5y (k) 15 + + 22 22 11

(1.5)

se sustituirá en las ecuaciones (1.5) las siguientes condiciones iniciales.

h

x (0) = 0

y (0) = 0

z (0) = 0

i

para obtener

h

x (1) =?

y (1) =?

z (1) =?

i

El proceso se realizara consecutivamente hasta que la norma entre dos vectores consecutivos sea menor a cierta tolerancia preestablecida.

Auxiliar: Martin Laguna A.

7 1.2. GAUSS SEIDEL

La norma θ se calcula con la ecuación: θ=

q

¡

x (k+1) − x (k)

¢2

¡ ¢2 ¡ ¢2 + y (k+1) − y (k) + z (k+1) − z (k) + ....+ ⇒ θ < tolerancia

(1.6)

Para proceder en la solución se tomara una tolerancia (10−5 ) , de esta forma realizaremos algunos cálculos iterativos con las ecuaciones (1.5) y (1.6), para luego mostrar los resultados en una tabla resumen de las iteraciones realizadas. Primera Iteración: h

Para k = 0 determinamos

x (1)

y (1)

x (1)

z (1)

=

=

29,411765

5x (0) 2z (0) 200 5 · 0 2 · 0 200 + + = + + 21 21 21 21 21 21

9,52381

=

z (1)

:

2y (0) 3z (0) 500 2 · 0 3 · 0 500 + + = + + 17 17 17 17 17 17

=

y (1)

i

=

5x (0) 5y (0) 15 5 · 0 5 · 0 15 + + = + + 22 22 11 22 22 11

=

1,363636

por tanto la norma θ1 :

θ1

q =

q =

=

¡

x (1) − x (0)

¢2

¡ ¢2 ¡ ¢2 + y (1) − y (0) + z (1) − z (0)

(29,411765 − 0)2 + (9,52381 − 0)2 + (1,363636 − 0)2

30,945345

θ1 < tolerancia ⇒ 30,945345 < 0,00001

Segunda Iteración: Para k = 1 determinamos

h

x (2)

y (2)

z (2)

i

:

Auxiliar: Martin Laguna A.

¡

¢

No se cumple

8 1.2. GAUSS SEIDEL

x (2)

=

2y (1) 3z (1) 500 2 · 9,52381 3 · 1,363636 500 + + = + + 17 17 17 17 17 17

30,772855

=

y (2)

5x (1) 2z (1) 200 5 · 29,411765 2 · 1,363636 200 + + = + + 21 21 21 21 21 21

=

16,656481

=

z (2)

=

5x (1) 5y (1) 15 5 · 29,411765 5 · 9,52381 15 + + = + + 22 22 11 22 22 11

=

10,212631

por tanto la norma θ2 :

θ2

q =

q =

=

¡

x (2) − x (1)

¢2

¡ ¢2 ¡ ¢2 + y (2) − y (1) + z (2) − z (1)

(30,772855 − 29,411765)2 + (16,656481 − 9,52381)2 + (10,212631 − 1,363636)2

11,446932

θ2 < tolerancia ⇒ 11,446932 < 0,00001

N

¡

¢

No se cumple

Repetimos el proceso hasta que la norma sea menor a la tolerancia aplicada, se presenta en la siguiente tabla todas las iteraciones realizadas.

Auxiliar: Martin Laguna A.

9 1.2. GAUSS SEIDEL

i

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

xi

yi

zi

θi

θi < tol

29.411765 30.772855 33.17358 33.651513 33.883473 33.955412 33.982015 33.991278 33.994536 33.995689 33.996094 33.996237 33.996287 33.996305 33.996311 33.996313

9.52381 16.656481 17.823311 18.57876 18.769773 18.851699 18.877982 18.887647 18.890997 18.892182 18.8926 18.892747 18.892799 18.892817 18.892823 18.892826

1.363636 10.212631 12.143031 12.953839 13.234153 13.330283 13.365252 13.377272 13.381574 13.383076 13.383607 13.383794 13.38386 13.383883 13.383891 13.383894

30.945345 11.446932 3.294149 1.206869 0.410935 0.145355 0.051199 0.017991

No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple No se cumple Se cumple

6,351696 × 10−3 2,234043 × 10−3 7,877488 × 10−4 2,77362 × 10−4 9,773312 × 10−5 3,442529 × 10−5 1,212773 × 10−5 4,272278 × 10−6

Tabla. 1.3: 16 iteraciones por el método Gauss-Seidel

Solucion: La solución del sistema por métodos iterativos: x (16) = 33,996 ; y (16) = 18,893 ; z (16) = 13,384

Una manera de verificar la solución es remplazar los resultados en el sistema matricial. Ecuación N° 1: 17x (16) − 2y (16) − 3z (16) = 500 17 · 33,996313 − 2 · 18,892826 − 3 · 13,383894 = 500 499,99999 = 500

Ecuación N° 2: −5x (16) + 21y (16) − 2z (16) = 200 −5 · 33,996313 + 21 · 18,892826 − 2 · 13,383894 = 200 199,99999 = 200

Ecuación N° 3: −5x (16) − 5y (16) + 22z (16) = 30 −5 · 33,996313 − 5 · 18,892826 + 22 · 13,383894 = 30 29,99997 = 30

Auxiliar: Martin Laguna A.

10 1.3. ELEMENTOS FINITOS EN UNA DIMENSIÓN

1.3

ELEMENTOS FINITOS EN UNA DIMENSIÓN

Ejercicio 3 Utilizando dos elementos finitos lineales calcular los desplazamientos, giros y dibujar el diagrama de desplazamiento. L=2 m

;

F = 18 kN

;

q = 10 kN

;

E = 210 GP a

;

I = 4 × 10−4 m 4

Desarrollo y solución

C

omo paso inicial realizaremos la malla del sistema estructural, como se observa en la figura 1.3.

fig. 1.3: Discretización de los elementos

Con ayuda de la fig. 1.3 determinamos el vector desplazamiento y el vector fuerza del sistema, para luego incorporar los mismo en el sistema mecánico estructural general. Auxiliar: Martin Laguna A.

11 1.3. ELEMENTOS FINITOS EN UNA DIMENSIÓN

1 1 2

n→ − o D est =

3 4 5



v1

         

 θ2   θ3    0   0   0

6



(1.7)

1 1 2

n→ − o F est =

3 4 5

          

6

0



 0   0    R4   R5   M6

(1.8)

vector desplazamiento-vector fuerza

Ahora pasamos a determinar la matriz de rigidez de cada miembro. miembro 1: 2E I = 168000 kN − m 2 , L 1 = 4 [m] £

¤



5

[k 1 ] =

5

6

1

2

31500

63000

−31500

63000

168000

−63000

  63000   −31500  63000

6 1 2

−63000

31500

84000

−63000



 84000   −63000   168000

miembro 2: E I = 84000 kN − m 2 , L 2 = 2 [m] £

¤

1

[k 2 ] =

2 4 3



1

2

4

3

126000

126000

−126000

126000

168000

−126000

−126000

126000

84000

−126000

  126000   −126000  126000



   −126000   168000 84000

Escribimos la matriz de conectividad [B ] que sera de mucha ayuda para ensamblar las matrices [ki ].

elemento 1 elemento 2

|

"

nodo 1

nodo 1

nodo 2

nodo 2

5

6

1

2

1

2 {z

4

3

Matriz −→ [B ]

Auxiliar: Martin Laguna A.

#

}

12 1.3. ELEMENTOS FINITOS EN UNA DIMENSIÓN

1 2

[K est ] =

3 4 5 6



1

2

3

4

5

6

157500

63000

126000

−126000

−31500

−63000

336000

84000

−126000

63000

84000

168000

−126000

0

−126000

−126000

126000

0

63000

0

0

31500

84000

0

0

63000

  63000   126000    −126000   −31500  −63000



 84000    0    0  63000   168000

matriz de rigidez de toda la estructura

Estudiamos la solicitación de las cargas en la estructura, por separado en cada elemento.

Elemento 1: carga concentrada

fig. 1.4: Descripción

de la fig. 1.4 determinamos el vector de carga: 1 5

© ª f1 =

6 1 2



−9



   −9     −9    9

Elemento 2: carga distribuida

fig. 1.5: Descripción

Auxiliar: Martin Laguna A.

(1.9)

13 1.3. ELEMENTOS FINITOS EN UNA DIMENSIÓN

de la fig. 1.5 determinamos el vector de carga: 1 1

© ª f2 =

2 4 3





−10

   −3,333     −10    3,333 → −

A estos dos vectores los ensamblamos para obtener Q est . 1 1 2

h→ − i Q est =

3 4 5 6



−19



   5,6667     3,3333       −10     −9    −9

(1.10)

vector de carga de toda la estructura

Con las ecuaciones (1.7),(1.8),(1.9) y (1.10) ya estamos en condiciones de plantear el sistema matricial estructural.

1 2 3 4 5 6



1

2

3

4

5

6

157500

63000

126000

−126000

−31500

−63000

336000

84000

−126000

63000

84000

168000

−126000

0

−126000

−126000

126000

0

63000

0

0

31500

84000

0

0

63000

  63000   126000    −126000   −31500  −63000

1



1

 84000   2  3 0  ·  4 0  63000   5 168000 6

1



v1

         

 θ2   θ3   = 0   0   0



1 2 3 4 5

          

6

0

1



1

 0   2 0   3 + R4   4 R5   5 M6 6



−19



   5,6667     3,3333       −10     −9    −9

(1.11) El sistema se reduce a un sistema matricial de 3 × 3, el cual se denota con los términos subrayados a continuación.             

 





 



  63000 126000 −126000 −31500 −63000  0  −19  157500  v1                                      θ 63000 336000 84000 −126000 63000 84000   0 5,6667       2                          84000 168000 −126000 0 0 0 3,3333  126000 θ3  = + ·          −126000 −126000 −126000 126000 0 0 R4  −10  0                                    −31500 63000 0 0 31500 63000   R −9      0 5                        −63000 84000 0 0 63000 168000 M6 −9  0

solución general del sistema (1.11) se puede obtener de manera directa utilizando un programa académico como lo es Mathcad que es de libre descarga, para simplificar el problema a continuación mostramos los resultados obtenidos de forma directa.

Auxiliar: Martin Laguna A.

14 1.3. ELEMENTOS FINITOS EN UNA DIMENSIÓN

Solucion:

Desplazamiento lineal: v 1 = −0,00034127 [m] v 4 = 0 [m] v 5 = 0 [m]

Desplazamiento rotacional: θ3 = 0,00026899 [r ad ] θ6 = 0 [r ad ] θ2 = 0,00001361 [r ad ]

N

Para realizar el diagrama de desplazamiento se tiene una solución cúbica supuesta. Dado que tanto el desplazamiento como la rotación se usan como grados de libertad para cada nodo, la solución supuesta apropiada se escribe utilizando la interpolación de Hermite. Las expresiones más convenientes se obtienen al introducir el siguiente cambio en las coordenadas:

s = x − x1 µ v (s) =

x1 ≤ x ≤ x2

2s 3 3s 2 − 2 +1 L3 L

0≤s ≤L s 3 3s 2 − +s L2 L

3s 2 2s 3 − 3 L2 L

  v1    ¶  θ 3 2 s s 1 · −  v L2 L   2   θ 2

          

v (s) = N T · d

(1.12)

Solución: para el elemento 1 en 0 ≤ x ≤ 4

Nodo inicial:

x1 = 0

Coordenadas locales:

s=x

Nodo final: Longitud:

x2 = 4 L=4

Usando estos valores y las funciones de interpolación en terminos de x se tiene: µ 3 ¶ x 3x 2 x3 x2 3x 2 x 3 x 3 x 2 NT = − +1 − +x − − 16 16 2 16 32 16 4´ ³ 32 T d = 0 0 −0,0003412698413 0,00001360544218

Auxiliar: Martin Laguna A.

(1.13)

15 1.3. ELEMENTOS FINITOS EN UNA DIMENSIÓN

remplazando los cálculos obtenidos en la ecuación (1.13) en la ecuación (1.12) se tiene la primera función: v 1 (x) = −0,000067389455789 · x 2 + 0,000011515022677 · x 3

(1.14)

Solución: para el elemento 2 en 4 ≤ x ≤ 6 De la misma forma se trabaja con el segundo elemento y se obtiene la siguiente función. v 2 (x) = 0,00226870748382 − 0,00155328798242 · x + 0,000283871882205 · x 2 + − 0,000014668367355 · x 3

(1.15)

Finalmente con la ecuaciones (1.14) y (1.15) graficamos el diagrama de desplazamiento de la estructura, para esto utilizamos el programa Mathcad como se ve en la figura 1.6.

fig. 1.6: Diagrama de desplazamiento aproximado

Para obtener resultados mas precisos en el gráfico utilizaremos un programa académico conocido como Ftool.

fig. 1.7: Diagrama de desplazamiento vertical

Auxiliar: Martin Laguna A.

16 1.3. ELEMENTOS FINITOS EN UNA DIMENSIÓN

fig. 1.8: Diagrama de desplazamiento rotacional

Auxiliar: Martin Laguna A.

17

2 2.1

Segundo Parcial

SERIES DE FOURIER

Ejercicio 1 Desarrolle la siguiente función en series de Fourier:    3 f (x) = 2x   1

;

0