Initia

Ejemplo de algoritmo de actualización Curso de Geomecánica Computacional Dr. Alejo O. Sfriso

Introducción Este ejercicio presenta el concepto básico de una implementación robusta de un algoritmo numérico. Se trata de encontrar el valor de y que cumple y = f[y] x con x dado. Por ejemplo In[19]:=

Clear  f , a , m , y  f [ y_ ] := a * ym

Solución analítica 1

El problema tiene dos soluciones analíticas: y = (a x) 1-m e y = 0. Por ejemplo In[21]:=

a = 1000 ; m = 0.5 ; Print  ( a * x0 )

1 1-m

x0 = 0.1 ; 

10 000.

Solución numérica, Mathematica Mathematica tiene incorporado el cálculo numérico de la solución, que depende del valor semilla que se elija como punto de partida. Por ejemplo, si se parte de 5000 o 1000 se llega respectivamente a In[23]:=

FindRoot  y - f [ y ] * x0 ⩵ 0 , { y , 5000 }  FindRoot  y - f [ y ] * x0 ⩵ 0 , { y , 1000 } 

Out[23]=

{y → 10 000.}

Out[24]=

y → 1.1659 × 10-23 

donde ambas soluciones son correctas, aunque una de las dos no interesa. Si se impone un cambio de x a xnew , el nuevo valor de y es In[29]:=

Clear [ y ] Δx = 0.02 ; xnew = x0 + Δx ; FindRoot  y - f [ y ] * xnew ⩵ 0 , { y , 50000 } 

2

2 Newton-Raphson.nb

Out[32]=

{y → 14 400.}

Que coincide con la solución analítica 1

In[34]:=

Print  ( a * xnew ) 1-m  14 400.

Solución numérica explícita Una manera tentadora - pero incorrecta - de resolver el problema consiste en convertir la expresión y = f[y] x en una forma incremental y = f[y] x y aplicar la fórmula yn+1 = yn + f[yn ] Δx que es un algoritmo explícito en x. En el código que sigue se demuestra su funcionamiento. y0 = 10000 ; Do y = y0 ; n = 10 ^ j ; dx = Δx / n ; Do  y = y + f [ y ] * dx ,  i , 1 , n   ; Print  "Cant. pasos=" , n , " Valor final =" , y  ; ,  j , 0 , 5   Cant. pasos=1 Valor final =12 000. Cant. pasos=10 Valor final =12 089.6 Cant. pasos=100 Valor final =12 099. Cant. pasos=1000 Valor final =12 099.9 Cant. pasos=10 000 Valor final =12 100. Cant. pasos=100 000 Valor final =12 100.

Como puede verse, el algoritmo no llega al valor correcto para ningún número de pasos.

Algoritmo de Newton-Raphson El problema puede resolverse correctamente con un esquema de Newton - Raphson convencional. El código es In[63]:=

Clear  x , y , Residuo , DerivadaResiduo  Residuo [ x_ , y_ ] := y - f [ y ] * x DerivadaResiduo [ x_ , y_ ] = D  Residuo [ x , y ] , y  ;

2 Newton-Raphson.nb

In[66]:=

3

Clear [ NewtonRaphson ] NewtonRaphson [ x_ , y0_ ] := Module   i = 1 , dy , y = y0 , output = Table  _ ,  i , 10    , output [ [ 1 ] ] = y0 ; While  Abs  Residuo [ x , y ]  > 0.0001 ∧ i < 10 , dy = Re  - Residuo [ x , y ]  DerivadaResiduo [ x , y ]  ; y = y + dy ; output   i + 1   = y ; i += 1 ;  ; Print [

output ] ;



La cantidad de pasos necesarias para alcanzar convergencia depende del valor semilla. Por ejemplo In[68]:=

NewtonRaphson [ xnew , 14000 ] NewtonRaphson [ xnew , 10000 ] NewtonRaphson [ xnew , 4000 ] {14 000, 14 402.9, 14 400., _, _, _, _, _, _, _} {10 000, 15 000., 14 405.9, 14 400., 14 400., _, _, _, _, _} {4000, 73 947.3, 20 935.1, 14 831.9, 14 403.1, 14 400., _, _, _, _}

El mismo algoritmo falla cuando el valor semilla está fuera del rango de convergencia del método. In[62]:=

NewtonRaphson [ xnew , 1000 ] {1000, - 2114.37, 1332.03, - 3400.51, 1748.71, - 5770.53, 2216.94, - 10 298.9, 2667.55, - 19 164.2}

Newton - Raphson + Bisección El algoritmo de Newton - Raphson se puede complementar con una bisección para obtener convergencia segura y rápida en el campo real

4

2 Newton-Raphson.nb

In[72]:=

Clear  NewtonMejorado  NewtonMejorado [ x_ , y0_ ] := Module   i = 1 , y = y0 , ymin = 10. ^ -16 , ymax = 10. ^ 16 , dy , yNR  , Print  "iter=0 {ymin,y,ymax}=" ,  ymin , y0 , ymax   ; While  Abs  Residuo [ x , y ]  > 0.01 ∧ i < 10 , (* Calcula el incremento de y según Newton-Raphson *) dy = Re  - Residuo [ x , y ]   Re  DerivadaResiduo [ x , y ]  ; (* Calcula y actualizado y lo guarda en una variable auxiliar yNR *) yNR = y + dy ; (*

Si yNR está en rango, lo acepta. Si está fuera de rango, lo descarta y elige un punto dentro del rango *)

If  yNR > ymin ∧ yNR < ymax , y = yNR , y = Sqrt  ymin * ymax   ; (* Evalúa la función y ajusta uno de los dos bordes de rango *) If  Residuo [ x , y ] > 0 , ymax = y , ymin = y  ; Print  "iter=" , i , " {ymin,y,ymax}=" ,  ymin , y , ymax   ; i += 1 ;  ; 

In[74]:=

NewtonMejorado [ xnew , 1000 ] iter=0 {ymin,y,ymax}=1. × 10-16 , 1000, 1. × 1016  iter=1 {ymin,y,ymax}=1., 1., 1. × 1016  iter=2 {ymin,y,ymax}=1., 1. × 108 , 1. × 108  iter=3 {ymin,y,ymax}={1., 603 622., 603 622.} iter=4 {ymin,y,ymax}={1., 50 517.1, 50 517.1} iter=5 {ymin,y,ymax}={1., 18 396.6, 18 396.6} iter=6 {ymin,y,ymax}={1., 14 593.9, 14 593.9} iter=7 {ymin,y,ymax}={1., 14 400.6, 14 400.6} iter=8 {ymin,y,ymax}={1., 14 400., 14 400.}