I implemented the derivatives for sec and cdc at https://github.com/sympy/sympy/pull/1716.
Aaron Meurer On Fri, Jan 18, 2013 at 4:12 PM, Aaron Meurer <[email protected]> wrote: > (hint, you can debug nsolve() by adding verbose=True) > > It looks like the issue is that sec doesn't have the derivative defined. > > Unfortunately, I tried defining it, but I then got > > ZeroDivisionError: matrix is numerically singular > > So there may be some more work for you to do. > > Aaron Meurer > > On Fri, Jan 18, 2013 at 1:30 PM, Renxian <[email protected]> wrote: >> why would the code bellow yield this error?????? thanks >> >> from sympy import * >> import numpy as np >> init_printing(use_unicode=False, wrap_line=True, no_global=True) >> #flat earth, rigid ,symatric , body axis >> theta,phi,psi=var('theta,phi,psi') >> p,q,r=var('p,q,r') >> u,v,w=var('u,v,w') >> Xe,Ye,Ze=var('Xe,Ye,Ze') >> deltAileron,deltRudder=var('deltAileron,deltRudder') >> deltFlap,deltElevator,dotAlpha=var('deltFlap,deltElevator,dotAlpha') >> >> #################################### >> #g=var('g') >> g=9.81#m/s**2 >> >> H_I2B=Matrix([[cos(theta)*cos(psi),cos(theta)*sin(psi),-sin(theta)], >> >> [(-cos(theta)*sin(psi)+sin(phi)*sin(theta)*cos(psi)),(cos(phi)*cos(psi)+sin(phi)*sin(theta)*cos(psi)),sin(phi)*cos(theta)], >> >> [(sin(phi)*sin(psi)+cos(phi)*sin(theta)*cos(psi)),(-sin(phi)*cos(psi)+cos(phi)*sin(theta)*sin(psi)),cos(phi)*cos(theta)]]) >> >> gx,gy,gz=H_I2B*Matrix(3,1,[0,0,g])#gravity in body frame >> >> #b,c,S=var('b,c,S') >> #wing span >> b=2.8956# m >> #mean aerodynamics chord >> c=0.189941 # % >> #wing area >> S=0.55# m**2 >> #m=var('m') >> m=13.5# kg >> #e,AR=var('e,AR') >> #Oswald's coefficient >> e=0.75 >> AR=b**2/S >> >> #CL_0,CL_alpha,CL_deltFlap,CL_deltElevator,CL_dotAlpha,CL_q,CL_Mach=var('CL_0,CL_alpha,CL_deltFlap,CL_deltElevator,CL_dotAlpha,CL_q,CL_Mach') >> ####Lift coefficient###### >> #Zero-alpha lift >> CL_0=0.23 >> CL_alpha=5.6106 >> CL_deltFlap=0.74 >> CL_deltElevator=0.13 >> CL_dotAlpha=1.9724 >> CL_q=7.9543 >> CL_Mach=0 >> #CD_0,CL_minD,CD_deltFlap,CD_deltElevator,CD_deltAileron,CD_deltRudder,CD_Mach=var('CD_0,CL_minD,CD_deltFlap,CD_deltElevator,CD_deltAileron,CD_deltRudder,CD_Mach') >> CD_0=0.0434 >> CL_minD=0.0434 >> CD_deltFlap=0.1467 >> CD_deltElevator=0.0135 >> CD_deltAileron=0.0302 >> CD_deltRudder=0.0303 >> CD_Mach=0 >> #CY_beta,CY_deltAileron,CY_deltRudder,CY_p,CY_r=var('CY_beta,CY_deltAileron,CY_deltRudder,CY_p,CY_r') >> CY_beta=-0.83 >> CY_deltAileron=-0.075 >> CY_deltRudder=0.1914 >> CY_p=0 >> CY_r=0 >> #Cm_0,Cm_alpha,Cm_deltFlap,Cm_deltElevator,Cm_dotAlpha,Cm_q,Cm_Mach=var('Cm_0,Cm_alpha,Cm_deltFlap,Cm_deltElevator,Cm_dotAlpha,Cm_q,Cm_Mach') >> Cm_0=0.135 >> Cm_alpha=-2.7397 >> Cm_deltFlap=0.0467 >> Cm_deltElevator=-0.9918 >> Cm_dotAlpha=-10.3796 >> Cm_q=-38.2067 >> Cm_Mach=0 >> #Cl_beta,Cl_deltAileron,Cl_deltRudder,Cl_p,Cl_r=var('Cl_beta,Cl_deltAileron,Cl_deltRudder,Cl_p,Cl_r') >> Cl_beta=-0.13 >> Cl_deltAileron=-0.1695 >> Cl_deltRudder=0.0024 >> Cl_p=-0.5051 >> Cl_r=0.2519 >> #Cn_beta,Cn_deltAileron,Cn_deltRudder,Cn_p,Cn_r=var('Cn_beta,Cn_deltAileron,Cn_deltRudder,Cn_p,Cn_r') >> Cn_beta=0.0726 >> Cn_deltAileron=0.0108 >> Cn_deltRudder=0.-0.0693 >> Cn_p=-0.069 >> Cn_r=-0.0946 >> >> ################ISA for Troposhere ############# >> T0=288.15 # sea level tempreture ,k (15 centigree) >> P_atm0=101325# N/m**2 >> R=287.04 #Characteristic gas constant (J/Kg/K) >> gamma=1.4 #Ratio of specific heats >> h=1000 >> T=T0-6.5*h/1000 >> P_atm=P_atm0*(1-0.0065*h/T0)**5.2561 >> rho=P_atm/(R*T) >> speedofsound=(T*R*gamma)**0.5 >> ####################################33 >> >> #uw,vw,ww=var('uw,vw,ww')#wind speed >> uw=0 >> vw=0 >> ww=0 >> >> ua,va,wa=Matrix(3,1,[u,v,w])-H_I2B*Matrix(3,1,[uw,vw,ww]) >> #alpha,beta,Va,Mach=var('alpha,beta,Va,Mach') >> Va=(ua**2+va**2+wa**2)**0.5 >> alpha=atan(wa/ua) >> beta=asin(va/ua) >> Mach=Va/speedofsound >> q_bar=rho*Va**2/2 >> >> #CL,CD,CY_w,Cm,Cl,Cn=var('CL,CD,CY,Cm,Cl,Cn') >> CL=CL_0+CL_alpha*alpha+CL_deltFlap*deltFlap+CL_deltElevator*deltElevator+(CL_dotAlpha*dotAlpha+CL_q*q)*c/(2*Va)+CL_Mach*Mach >> #CD=CD_0+(CL-CL_minD)**2/(pi*e*AR)+CD_deltFlap*Abs(deltFlap)+CD_deltElevator*Abs(deltElevator)+CD_deltAileron*Abs(deltAileron)+CD_deltRudder*Abs(deltRudder)+CD_Mach*Mach >> CD=CD_0+(CL-CL_minD)**2/(pi*e*AR)+CD_deltFlap*(deltFlap)+CD_deltElevator*(deltElevator)+CD_deltAileron*(deltAileron)+CD_deltRudder*(deltRudder)+CD_Mach*Mach >> CY_w=CY_beta*beta+CY_deltAileron*deltAileron+CY_deltRudder*deltRudder+(CY_p*p+CY_r*r)*b/(2*Va) >> Cm=Cm_0+Cm_alpha*alpha+Cm_deltFlap*deltFlap+Cm_deltElevator*deltElevator+(Cm_dotAlpha*dotAlpha+Cm_q*q)*c/(2*Va)+Cm_Mach*Mach >> Cl=Cl_beta*beta+Cl_deltAileron*deltAileron+Cl_deltRudder*deltRudder+(Cl_p*p+Cl_r*r)*b/(2*Va) >> Cn=Cn_beta*beta+Cn_deltAileron*deltAileron+Cn_deltRudder*deltRudder+(Cn_p*p+Cn_r*r)*b/(2*Va) >> H_b2w=Matrix([[cos(alpha)*cos(beta),sin(beta),sin(alpha)*cos(beta)], >> [-cos(alpha)*sin(beta),cos(beta),-sin(alpha)*sin(beta)], >> [-sin(alpha),0,cos(alpha)]]) >> #CX,CY,CZ=var('CX,CY,CZ') >> CX,CY,CZ=H_b2w*Matrix(3,1,[-CD,CY_w,-CL]) >> >> >> ##############Fixed-pitch propeller##############333 >> >> #Jar,CT,CP=var('Jar,CT,CP') >> omega=var('omega') >> #propeller radius >> R_prop=0.254# m >> >> Jar=pi*Va/(omega*R_prop) >> #coefficient of thrust (CT) and power (CP) >> #the fomular is got after curve fitting using quadratic polynomial >> CT=-0.02162763*Jar**2-0.03491274*Jar+0.03787723 >> CP=-0.01860428*Jar**2-0.01512463*Jar+0.02509024 >> F_prop=4*rho*R_prop**4*omega**2*CT/pi**2 >> M_prop=-4*rho*R_prop**5*omega**2*CP >> ##################Piston engine################### >> thr=var('thr') >> MAPmin=60; >> MAP=thr*(P_atm/1000-MAPmin)+MAPmin >> RPM=omega*30/pi >> #the fomulas of fuel flow and power of piston are got using polynomial >> surface fitiing in Matlab >> Fuelflow=-3189+0.2615*RPM+119.6*MAP-0.00001329*RPM**2-0.005922*RPM*MAP-1.843*MAP**2+3.257/10**10*RPM**3+1.377/10**7*RPM**2*MAP+0.00003562*RPM*MAP**2+0.006149*MAP**3 >> Power_P=-10600+1.38*RPM+367.1*MAP-0.0000768*RPM**2-0.03094*RPM*MAP-4.202*MAP**2+3.727/10**9*RPM**3+3.917/10**7*RPM**2*MAP+0.0002021*RPM*MAP**2+0.01572*MAP**3 >> #Power_P=432.4-0.3574*RPM-4.54*MAP+0.000005236*RPM**2+0.004861*RPM*MAP >> #Fuelflow,Power_P=var('Fuelflow,Power_P') >> #Power_P=var('Power_P') >> power_corrected=(T0/T)**0.5*Power_P >> M_eng=power_corrected/omega >> ###### engine total force and moment######### >> J_eng=0.001# engine shaft moment of inertia >> J_prop=0.002# propeller moment of inertia >> MomentofInertiaofEngAndProp=(M_eng+M_prop)/(J_eng+J_prop) >> #omega=integrate(MomentofInertiaofEngAndProp,omega) >> #omega2=integrate(MomentofInertiaofEngAndProp,omega) >> #without consider the moment caused by proleller force >> Tx=F_prop >> Ty=0 >> Tz=0 >> M_Tx=-M_eng >> M_Ty=0 >> M_Tz=0 >> #############equations of motion################ >> #Ixx,Iyy,Izz,Ixz=var('Ixx,Iyy,Izz,Ixz') >> #gross moment of inertia >> Ixx=0.8244# kg*m**2 >> Iyy=1.135# kg*m**2 >> Izz=1.759# kg*m**2 >> Ixz=0.1204# kg*m**2 >> >> X=CX*q_bar*S+Tx >> Y=CY*q_bar*S+Ty >> Z=CZ*q_bar*S+Tz >> L=Cl*q_bar*S*b+M_Tx >> M=Cm*q_bar*S*c+M_Ty >> N=Cn*q_bar*S*b+M_Tz >> >> u_dot=X/m+gx+r*v-q*w >> v_dot=Y/m+gy-r*u+p*w >> w_dot=Z/m+gz+q*u-p*v >> p_dot=(Izz*L+Ixz*N-(Ixz*(Iyy-Ixx-Izz)*p+(Ixx**2+Izz*(Izz-Iyy))*r)*q)/(Ixx*Izz-Izz**2) >> q_dot=(M-(Ixx-Izz)*p*r-Ixz*(p**2-r**2))/Iyy >> r_dot=(Ixz*L+Ixx*N+(Ixz*(Iyy-Ixx-Izz)*r+(Ixz**2+Ixx*(Ixx-Iyy)*p)*q))/(Ixx*Izz-Ixz**2) >> x_I_dot=cos(theta)*cos(psi)*u+(-cos(phi)*sin(psi)+sin(phi)*sin(theta)*cos(psi))*v+(sin(phi)*sin(psi)+cos(phi)*sin(theta)*cos(psi))*w >> y_I_dot=cos(theta)*sin(psi)*u+(cos(phi)*cos(psi)+sin(phi)*sin(theta)*sin(psi))*v+(-sin(phi)*cos(psi)+cos(phi)*sin(theta)*sin(psi))*w >> z_I_dot=-sin(theta)*u+sin(phi)*cos(theta)*v+cos(phi)*cos(theta)*w >> phi_dot=p+(q*sin(phi)+r*cos(phi))*tan(theta) >> theta_dot=q*cos(phi)-r*sin(phi) >> psi_dot=(q*sin(phi)+r*cos(phi))*sec(theta) >> >> ############# Trim for straight and level flight ############ >> ###the aircraft is trimmed at >> Va25=Va-25 >> #deltFlap=0 >> #deltAlpha=0 >> #straight level fight means >> #phi=0 >> #p=0 >> #r=0 >> >> ###Initial seed value for solver >> #q=0 >> #deltElevator=0 >> #deltAileron=0 >> #deltRudder=0 >> #thr=0.5 >> #u=25 >> #v=0 >> #w=0 >> #omega=5236 >> #theta=0 >> #psi=0 >> ###solve >> u_dot_=u_dot.subs({deltElevator:0,deltAileron:0,deltRudder:0,deltFlap:0,dotAlpha:0,omega:545.4,thr:0.07063}) >> v_dot_=v_dot.subs({deltElevator:0,deltAileron:0,deltRudder:0,deltFlap:0,dotAlpha:0,omega:545.4,thr:0.07063}) >> w_dot_=w_dot.subs({deltElevator:0,deltAileron:0,deltRudder:0,deltFlap:0,dotAlpha:0,omega:545.4,thr:0.07063}) >> p_dot_=p_dot.subs({deltElevator:0,deltAileron:0,deltRudder:0,deltFlap:0,dotAlpha:0,omega:545.4,thr:0.07063}) >> q_dot_=q_dot.subs({deltElevator:0,deltAileron:0,deltRudder:0,deltFlap:0,dotAlpha:0,omega:545.4,thr:0.07063}) >> r_dot_=r_dot.subs({deltElevator:0,deltAileron:0,deltRudder:0,deltFlap:0,dotAlpha:0,omega:545.4,thr:0.07063}) >> phi_dot_=phi_dot.subs({deltElevator:0,deltAileron:0,deltRudder:0,deltFlap:0,dotAlpha:0,omega:545.4,thr:0.07063}) >> theta_dot_=theta_dot.subs({deltElevator:0,deltAileron:0,deltRudder:0,deltFlap:0,dotAlpha:0,omega:545.4,thr:0.07063}) >> psi_dot_=psi_dot.subs({deltElevator:0,deltAileron:0,deltRudder:0,deltFlap:0,dotAlpha:0,omega:545.4,thr:0.07063}) >> #omega_=omega.subs([p,r,phi,deltFlap,deltAlpha],[0,0,0,0,0]) >> >> #result=nsolve([Va25,u_dot,v_dot,w_dot,p_dot,q_dot,r_dot,phi_dot,theta_dot,psi_dot,omega], >> # >> [q,deltElevator,deltAileron,deltRudder,thr,u,v,w,omega,theta,psi], >> # [0,0,0,0,0.5,25,0,0,5236,0,0]) >> #result=nsolve([u_dot_,v_dot_,w_dot_,p_dot_,q_dot_,r_dot_,phi_dot_,theta_dot_,psi_dot_], >> # [q,deltElevator,deltAileron,deltRudder,u,v,w,theta,psi], >> # >> [0.000,0.0000,0.00000,0.0000,25.0,0.0000,0.00000,0.00000,0.00000]) >> result=nsolve([u_dot_,v_dot_,w_dot_,p_dot_,q_dot_,r_dot_,phi_dot_,theta_dot_,psi_dot_], >> [q,r,p,phi,u,v,w,theta,psi], >> >> [0.000,0.0000,0.00000,0.0000,2.0,0.0000,0.00000,0.00000,0.00000]) >> >> >> -- >> You received this message because you are subscribed to the Google Groups >> "sympy" group. >> To view this discussion on the web visit >> https://groups.google.com/d/msg/sympy/-/IEdPzm9B8dYJ. >> To post to this group, send email to [email protected]. >> To unsubscribe from this group, send email to >> [email protected]. >> For more options, visit this group at >> http://groups.google.com/group/sympy?hl=en. -- You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
