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 originalfrom pylab import *
#insertamos la funcion con las condiciones oportunas
g=9.81
l1=3
l2=7
m1=10000
m2=200
r=0.75
c=2 # medidas de m1
b=1
mb=50 #masa de la barra
def fcn (x, y):
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)])
#def errg (x, y):
#return array ([cos(x), -sin(x)]) - y
#def errorl (x, y, xn, yn):
# return array ([yn[0]*cos(x-xn)+yn[1]*sin(x-xn), -yn[0]*sin(x-xn)+yn[1]*cos(x-xn)])-y
#Definimos los coeficientes de la matriz RK
#c es un vector de dimension 7
c2=1.0E0/5.0E0
c3=3.0E0/10.0E0
c4=4.0E0/5.0E0
c5=8.0E0/9.0E0
c6=1.0E0
c7=1.0E0
#La matriz a es de dimension 7
a21=1.0E0/5.0E0
a31=3.0E0/40.0E0
a41=44.0E0/45.0E0
a51=19372.0E0/6561.0E0
a61=9017.0E0/3168.0E0
a71=35.0E0/384.0E0
a32=9.0E0/40.0E0
a42=-56.0E0/15.0E0
a52=-25360.0E0/2187.0E0
a62=-355.0E0/33.0E0
a43=32.0E0/9.0E0
a53=64448.0E0/6561.0E0
a63=46732.0E0/5247.0E0
a73=500.0E0/1113.0E0
a54=-212.0E0/729.0E0
a64=49.0E0/176.0E0
a74=125.0E0/192.0E0
a65=-5103.0E0/18656.0E0
a75=-2187.0E0/6784.0E0
a76=11.0E0/84.0E0
#el vector b es el del paso mayor (orden5)
b1=35.0E0/384.0E0
b3=500.0E0/1113.0E0
b4=125.0E0/192.0E0
b5=-2187.0E0/6784.0E0
b6=11.0E0/84.0E0
# El vector be es el de orden menor (orden 4)
be1=5179.0E0/57600.0E0
be2=0
be3=7571.0E0/16695.0E0
be4=393.0E0/640.0E0
be5=-92097.0E0/339200.0E0
be6=187.0E0/2100.0E0
be7=1.0E0/40.0E0
numrecha=0
numacept=0
#intervalo en el que hay que realizar los calculos
x0=0.0; xf=5.0
#errores abs y rel
EABS=1e-8;EREL=1e-9
y0=array([arccos(-1.0)/4, 0.0]);
y=zeros([size(y0)])
n=1000;
tol=EABS+EREL*norm(y0)
g1=fcn(x0,y0)
h=min(0.5, (tol/(1e-38+norm(g1)))**(1/5))
print h
xab=zeros([n,1]);
yab=zeros([n,2]);
zab=zeros([n,2]);
fallo=False
ydibujo=zeros([10000, 2])
while x0<xf:
g2=fcn(x0+c2*h,y0+h*a21*g1)
g3=fcn(x0+c3*h,y0+h*(a31*g1+a32*g2))
g4=fcn(x0+c4*h,y0+h*(a41*g1+a42*g2+a43*g3))
g5=fcn(x0+c5*h,y0+h*(a51*g1+a52*g2+a53*g3+a54*g4))
g6=fcn(x0+c6*h,y0+h*(a61*g1+a62*g2+a63*g3+a64*g4+a65*g5))
g7=fcn(x0+c7*h,y0+h*(a71*g1+a73*g3+a74*g4+a75*g5+a76*g6))
#Calculamos la y
y5=y0+h*(b1*g1+b3*g3+b4*g4+b5*g5+b6*g6)
y4=y0+h*(be1*g1+be2*g2+be3*g3+be4*g4+be5*g5+be6*g6+be7*g7)
#tolerancia y estimación
est=norm(y5-y4)
tol=EABS+EREL*max(1,norm(y5))
if est<= tol:
#AQUI EL PASO SE ACEPTA
x0=x0+h
y0=y5
ydibujo[numacept,:]=y5
g1=g7
fac=min(2,0.9*(tol/(est+1e-38))**(1/5))
if fallo:
fac=min(1,fac)
h=fac*h
if x0+h>xf:
h=xf-x0
fallo=False
#print " angul ", x0, y5
numacept=numacept+1
else:
#El paso ha sido rechazado por tanto hay que recalcular el paso.
h=max(0.1, 0.9*(tol/est)**(1/5))*h
if fallo:
print "Dos pasos rechazados consecutivos"
fallo=True
numrecha=numrecha+1
print numrecha
print numacept
plot()
#plot (ydibujo[0:numacept,0],'g-')
show()