d2jsp
Log InRegister
d2jsp Forums > Off-Topic > Computers & IT > Programming & Development > Help With Python > From Scipy.optimize Import Minimize
Add Reply New Topic New Poll
Member
Posts: 7,459
Joined: Jul 31 2009
Gold: 3,618.00
Nov 7 2017 02:47pm
Ok so I'm trying to minimize a function, but it seems like its not even trying to. Not to mention that when I print ki2 its always 0 no matter what. The .txt im opening contains values for t,x, siigx,y,sigy

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

oméga=np.pi/4
w=0
i=np.pi/4
a=1
Tabe=0
Me=1.2
Mo=0
PaPa=(a,Mo,Tabe,oméga,w,i)
t,x,sigx,y,sigy=np.loadtxt("mesuresOrbiteQ4.txt",unpack=True)

def fct(E,M,e):
return M-E+e*np.sin(E)
def fctDer(E,e): #dérivée de la fonction fct
return -1+e*np.cos(E)
def newton(fct,fctDer,M,epsilon=1.e-10,iterMax=20):

prec=1. #Estimé de la précision, initialise à >epsilon
if M%(2*np.pi)>=0 and M%(2*np.pi)<np.pi:
Eo=M+Tabe/2
if M%(2*np.pi)>=np.pi and M%(2*np.pi)<2*np.pi:
Eo=M-Tabe/2

k=0 #compteur d'itérations
while (prec>epsilon) and (k<iterMax): #Boucle Newton
dx=-fct(Eo,M,Tabe)/fctDer(Eo,Tabe) #incrément en x
Eo+=dx #nouvelle valeur de x0 (la racine)
prec=abs(dx) #estimé conservateur de la précision
k+=1

return Eo

def position(PaPa,t,Me):

a=PaPa[0]
Mo=PaPa[1]
Tabe=PaPa[2]
oméga=PaPa[3]
w=PaPa[4]
i=PaPa[5]
Me=1.2

Tabx=np.zeros(len(t))
Taby=np.zeros(len(t))
tx=np.zeros(len(Tabx))
ty=np.zeros(len(Taby))

for j in range(len(t)):
P=np.sqrt(a**3/Me)
M=Mo+2*np.pi*t[j]/P
Eo=newton(fct,fctDer,M)
#M=Eo+Tabe*np.sin(Eo)
r=a*(1-Tabe*np.cos(Eo))
thêta=2*np.arctan2(np.sqrt(1-Tabe)*np.cos(Eo/2),np.sqrt(1+Tabe)*np.sin(Eo/2))

Tabx[j]=r*np.cos(thêta)
Taby[j]=r*np.sin(thêta)
tx[j]=(np.cos(w)*np.cos(oméga)-np.cos(i)*np.sin(oméga)*np.sin(w))*Tabx[j]-(np.sin(w)*np.cos(oméga)+np.cos(i)*np.sin(oméga)*np.cos(w))*Taby[j]
ty[j]=(np.cos(w)*np.sin(oméga)+np.cos(i)*np.cos(oméga)*np.sin(w))*Tabx[j]-(np.sin(w)*np.sin(oméga)-np.cos(i)*np.cos(oméga)*np.cos(w))*Taby[j]




return x,y,Tabx,Taby

def ki2fct(position,t,Me,x,sigx,y,sigy):

ki2=np.sum(((x-tx)/sigx)**2+((y-ty)/sigy)**2)

return ki2

t,x,sigx,y,sigy=np.loadtxt("mesuresOrbiteQ4.txt",unpack=True)

ki2=ki2fct(position,t,Me,x,sigx,y,sigy)
bornes=((0,1),(0.1,10),(0,5),(0,2*np.pi),(0,5),(0,2*np.pi))
p0=(a,Mo,Tabe,oméga,w,i)
res=minimize(ki2fct,p0,args=(t,Me,x,sigx,y,sigy),method="SLSQP",bounds=bornes)
tx,ty,tabx,taby=position(res.x,t,Me)


plt.plot(tx,ty,linestyle='None', marker='$\\clubsuit$',color='k', markersize=10)
plt.plot(tx,ty)
plt.title('Orbite',fontdict={'size': 14})
plt.plot([0,0],[1,-1],'--k',[-1,1],[0,0],'--k')
Go Back To Programming & Development Topic List
Add Reply New Topic New Poll