Foros del Web » Programando para Internet » Python »

Gráfica en mi Runge-Kutta.

 Estas en el tema de Gráfica en mi Runge-Kutta. en el foro de Python en Foros del Web. Hola ! Estoy programando una ecuación de un movimiento que realiza el lanzamiento de un proyectil con una catapulta con Runge-Kutta. El problema es que ...
  #1 (permalink)  
Antiguo 27/04/2016, 10:42
 
Fecha de Ingreso: abril-2016
Mensajes: 1
Antigüedad: 8 años, 6 meses
Puntos: 0
Pregunta Gráfica en mi Runge-Kutta.
Hola ! Estoy programando una ecuación de un movimiento que realiza el lanzamiento de un proyectil con una catapulta con Runge-Kutta.

El problema es que , tengo ya metida la ecuación en el array y hecho el runge kutta pero no se como pintar la solución. No se que poner , para que me pinte la solución.

Les dejo el programa, agradecería mucho la ayuda.

Código Python: 
Ver original
  1. from pylab import *
  2.  
  3. #insertamos la funcion con las condiciones oportunas
  4. g=9.81
  5. l1=3
  6. l2=7
  7.  
  8. m1=10000 
  9. m2=200
  10. r=0.75
  11. c=2 # medidas de m1
  12. b=1
  13. mb=50 #masa de la barra
  14.  
  15.  
  16. def fcn (x, y):
  17.  
  18.     return array ([y[1], -g*(m1*l1-m2*l2 + mb*(l1+l1)/2)*cos(y[0])/(m1*l1**2+m2*l2**2+(2/5)*m2*r**2+m1/12*(b**2+c**2) + 1/12*mb*(l1+l2)**2)])
  19.  
  20. #def errg (x, y):
  21.     #return array ([cos(x), -sin(x)]) - y
  22.  
  23. #def errorl (x, y, xn, yn):
  24.    # return array ([yn[0]*cos(x-xn)+yn[1]*sin(x-xn), -yn[0]*sin(x-xn)+yn[1]*cos(x-xn)])-y
  25.  
  26. #Definimos los coeficientes de la matriz RK
  27. #c es un vector de dimension 7 
  28. c2=1.0E0/5.0E0
  29. c3=3.0E0/10.0E0
  30. c4=4.0E0/5.0E0
  31. c5=8.0E0/9.0E0
  32. c6=1.0E0
  33. c7=1.0E0
  34.  
  35. #La matriz a es de dimension 7
  36. a21=1.0E0/5.0E0
  37. a31=3.0E0/40.0E0
  38. a41=44.0E0/45.0E0
  39. a51=19372.0E0/6561.0E0
  40. a61=9017.0E0/3168.0E0
  41. a71=35.0E0/384.0E0
  42. a32=9.0E0/40.0E0
  43. a42=-56.0E0/15.0E0
  44. a52=-25360.0E0/2187.0E0
  45. a62=-355.0E0/33.0E0
  46. a43=32.0E0/9.0E0
  47. a53=64448.0E0/6561.0E0
  48. a63=46732.0E0/5247.0E0
  49. a73=500.0E0/1113.0E0
  50. a54=-212.0E0/729.0E0
  51. a64=49.0E0/176.0E0
  52. a74=125.0E0/192.0E0
  53. a65=-5103.0E0/18656.0E0
  54. a75=-2187.0E0/6784.0E0
  55. a76=11.0E0/84.0E0
  56.       
  57. #el vector b es el del paso mayor (orden5)
  58. b1=35.0E0/384.0E0
  59. b3=500.0E0/1113.0E0
  60. b4=125.0E0/192.0E0
  61. b5=-2187.0E0/6784.0E0
  62. b6=11.0E0/84.0E0
  63.  
  64. # El vector be es el de orden menor (orden 4)
  65. be1=5179.0E0/57600.0E0
  66. be2=0
  67. be3=7571.0E0/16695.0E0
  68. be4=393.0E0/640.0E0
  69. be5=-92097.0E0/339200.0E0
  70. be6=187.0E0/2100.0E0
  71. be7=1.0E0/40.0E0
  72.     
  73. numrecha=0
  74. numacept=0
  75.  
  76. #intervalo en el que hay que realizar los calculos
  77. x0=0.0; xf=5.0
  78. #errores abs y rel
  79. EABS=1e-8;EREL=1e-9
  80. y0=array([arccos(-1.0)/4, 0.0]); 
  81. y=zeros([size(y0)])
  82. n=1000; 
  83.     
  84. tol=EABS+EREL*norm(y0)  
  85. g1=fcn(x0,y0)
  86. h=min(0.5, (tol/(1e-38+norm(g1)))**(1/5))
  87. print h
  88. xab=zeros([n,1]);
  89. yab=zeros([n,2]);
  90. zab=zeros([n,2]);
  91.  
  92. fallo=False 
  93. ydibujo=zeros([10000, 2])
  94. while x0<xf:
  95.     
  96.     g2=fcn(x0+c2*h,y0+h*a21*g1)
  97.     g3=fcn(x0+c3*h,y0+h*(a31*g1+a32*g2))
  98.     g4=fcn(x0+c4*h,y0+h*(a41*g1+a42*g2+a43*g3))
  99.     g5=fcn(x0+c5*h,y0+h*(a51*g1+a52*g2+a53*g3+a54*g4))
  100.     g6=fcn(x0+c6*h,y0+h*(a61*g1+a62*g2+a63*g3+a64*g4+a65*g5))
  101.     g7=fcn(x0+c7*h,y0+h*(a71*g1+a73*g3+a74*g4+a75*g5+a76*g6))
  102.     
  103.     
  104.     #Calculamos la y 
  105.     y5=y0+h*(b1*g1+b3*g3+b4*g4+b5*g5+b6*g6)
  106.     y4=y0+h*(be1*g1+be2*g2+be3*g3+be4*g4+be5*g5+be6*g6+be7*g7)
  107.     
  108.     #tolerancia y estimación
  109.     est=norm(y5-y4)
  110.     tol=EABS+EREL*max(1,norm(y5))
  111.  
  112.     if est<= tol:
  113.     #AQUI EL PASO SE ACEPTA  
  114.         x0=x0+h
  115.         y0=y5
  116.         ydibujo[numacept,:]=y5
  117.         
  118.         g1=g7
  119.         
  120.         fac=min(2,0.9*(tol/(est+1e-38))**(1/5))
  121.        
  122.         
  123.         if fallo:
  124.             fac=min(1,fac)
  125.         h=fac*h
  126.         if x0+h>xf:
  127.             h=xf-x0
  128.         fallo=False
  129.         
  130.        
  131.  
  132.         #print " angul ", x0, y5
  133.         numacept=numacept+1
  134.        
  135.     else:
  136.         
  137.         #El paso ha sido rechazado por tanto hay que recalcular el paso.
  138.         
  139.         h=max(0.1, 0.9*(tol/est)**(1/5))*h
  140.  
  141.         if fallo:
  142.             print "Dos pasos rechazados consecutivos"
  143.         fallo=True 
  144.         numrecha=numrecha+1
  145.  
  146. print numrecha
  147. print numacept
  148.  
  149. plot()
  150. #plot (ydibujo[0:numacept,0],'g-')
  151.  
  152. show()
Última edición por razpeitia; 08/05/2016 a las 10:33
  #2 (permalink)  
Antiguo 08/05/2016, 11:38
Avatar de razpeitia
Moderador
  
Fecha de Ingreso: marzo-2005
Ubicación: Monterrey, México
Mensajes: 7.321
Antigüedad: 19 años, 8 meses
Puntos: 1360
Respuesta: Gráfica en mi Runge-Kutta.
Tu script no me funcionaba, numacept hacia que accederia a ydibujo fuera de su rango.
Así que le agregue una condición para que parara una vez que excediera ese limite.
No le estabas mandado nada a plot, así que obviamente no ponía nada.

Imprimí los valores de h, x0 y xf y bueno h tiende a cero muy rápido lo que en 10000 iteraciones no iba a llevar a x0 a superar a xf.

Código Python: 
Ver original
  1. # coding: utf-8
  2. from pylab import *
  3.  
  4. # Insertamos la funcion con las condiciones oportunas
  5. g = 9.81
  6. l1 = 3
  7. l2 = 7
  8.  
  9. m1 = 10000
  10. m2 = 200
  11. r = 0.75
  12. c = 2  # medidas de m1
  13. b = 1
  14. mb = 50  # masa de la barra
  15.  
  16.  
  17. def fcn(x, y):
  18.     return array ([y[1], -g*(m1*l1-m2*l2 + mb*(l1+l1)/2)*cos(y[0])/(m1*l1**2+m2*l2**2+(2/5)*m2*r**2+m1/12*(b**2+c**2) + 1/12*mb*(l1+l2)**2)])
  19.  
  20.  
  21. # Definimos los coeficientes de la matriz RK
  22. # c es un vector de dimension 7
  23. c2 = 1.0E0 / 5.0E0
  24. c3 = 3.0E0 / 10.0E0
  25. c4 = 4.0E0 / 5.0E0
  26. c5 = 8.0E0 / 9.0E0
  27. c6 = 1.0E0
  28. c7 = 1.0E0
  29.  
  30. # La matriz a es de dimension 7
  31. a21 = 1.0E0 / 5.0E0
  32. a31 = 3.0E0 / 40.0E0
  33. a41 = 44.0E0 / 45.0E0
  34. a51 = 19372.0E0 / 6561.0E0
  35. a61 = 9017.0E0 / 3168.0E0
  36. a71 = 35.0E0 / 384.0E0
  37. a32 = 9.0E0 / 40.0E0
  38. a42 = -56.0E0 / 15.0E0
  39. a52 = -25360.0E0 / 2187.0E0
  40. a62 = -355.0E0 / 33.0E0
  41. a43 = 32.0E0/9.0E0
  42. a53 = 64448.0E0/6561.0E0
  43. a63 = 46732.0E0/5247.0E0
  44. a73 = 500.0E0/1113.0E0
  45. a54 = -212.0E0/729.0E0
  46. a64 = 49.0E0/176.0E0
  47. a74 = 125.0E0/192.0E0
  48. a65 = -5103.0E0/18656.0E0
  49. a75 = -2187.0E0/6784.0E0
  50. a76 = 11.0E0/84.0E0
  51.  
  52. # el vector b es el del paso mayor (orden5)
  53. b1 = 35.0E0 / 384.0E0
  54. b3 = 500.0E0 / 1113.0E0
  55. b4 = 125.0E0 / 192.0E0
  56. b5 = -2187.0E0 / 6784.0E0
  57. b6 = 11.0E0 / 84.0E0
  58.  
  59. # El vector be es el de orden menor (orden 4)
  60. be1 = 5179.0E0 / 57600.0E0
  61. be2 = 0
  62. be3 = 7571.0E0 / 16695.0E0
  63. be4 = 393.0E0 / 640.0E0
  64. be5 = -92097.0E0 / 339200.0E0
  65. be6 = 187.0E0 / 2100.0E0
  66. be7 = 1.0E0 / 40.0E0
  67.  
  68. numrecha = 0
  69. numacept = 0
  70.  
  71. # intervalo en el que hay que realizar los calculos
  72. x0 = 0.0
  73. xf = 5.0
  74.  
  75. # errores abs y rel
  76. EABS = 1e-8
  77. EREL = 1e-9
  78. y0 = array([arccos(-1.0)/4, 0.0])
  79. y = zeros([size(y0)])
  80. n = 1000
  81.  
  82. tol = EABS + EREL * norm(y0)
  83. g1 = fcn(x0, y0)
  84. h = min(0.5, (tol/(1e-38+norm(g1)))**(1/5))
  85. print h
  86. xab = zeros([n, 1])
  87. yab = zeros([n, 2])
  88. zab = zeros([n, 2])
  89.  
  90. fallo = False
  91. ydibujo = zeros([10000, 2])
  92. while x0 < xf and numacept < 10000:
  93.     g2 = fcn(x0+c2*h, y0+h*a21*g1)
  94.     g3 = fcn(x0+c3*h, y0+h*(a31*g1+a32*g2))
  95.     g4 = fcn(x0+c4*h, y0+h*(a41*g1+a42*g2+a43*g3))
  96.     g5 = fcn(x0+c5*h, y0+h*(a51*g1+a52*g2+a53*g3+a54*g4))
  97.     g6 = fcn(x0+c6*h, y0+h*(a61*g1+a62*g2+a63*g3+a64*g4+a65*g5))
  98.     g7 = fcn(x0+c7*h, y0+h*(a71*g1+a73*g3+a74*g4+a75*g5+a76*g6))
  99.  
  100.     # Calculamos la y
  101.     y5 = y0+h*(b1*g1+b3*g3+b4*g4+b5*g5+b6*g6)
  102.     y4 = y0+h*(be1*g1+be2*g2+be3*g3+be4*g4+be5*g5+be6*g6+be7*g7)
  103.  
  104.     # tolerancia y estimación
  105.     est = norm(y5-y4)
  106.     tol = EABS+EREL*max(1, norm(y5))
  107.  
  108.     if est <= tol:
  109.         # AQUI EL PASO SE ACEPTA
  110.         x0 = x0 + h
  111.         y0 = y5
  112.         ydibujo[numacept, :] = y5
  113.         g1 = g7
  114.         fac = min(2, 0.9*(tol/(est+1e-38))**(1/5))
  115.         if fallo:
  116.             fac = min(1, fac)
  117.         h = fac*h
  118.         if x0+h > xf:
  119.             h = xf-x0
  120.         fallo = False
  121.         numacept = numacept + 1
  122.     else:
  123.         # El paso ha sido rechazado por tanto hay que recalcular el paso.
  124.         h = max(0.1, 0.9*(tol/est)**(1/5))*h
  125.  
  126.         if fallo:
  127.             print "Dos pasos rechazados consecutivos"
  128.         fallo = True
  129.         numrecha = numrecha+1
  130.  
  131. print numrecha
  132. print numacept
  133.  
  134. plot(ydibujo)
  135. show()

También encontré una implementación usando funciones lambda.
https://rosettacode.org/wiki/Runge-Kutta_method#Python

Etiquetas: Ninguno
Atención: Estás leyendo un tema que no tiene actividad desde hace más de 6 MESES, te recomendamos abrir un Nuevo tema en lugar de responder al actual.
Respuesta


« Tema Anterior | Próximo Tema »


La zona horaria es GMT -6. Ahora son las 07:05.