0
|
9037 Views
|
|
0 Total Likes
View groups...
Share
GROUPS:

# How to solve algebraic looping in Mathematica ?

Posted 9 years ago
 Hi guys so here is the code (*Parameters*) tend = 100; Mu = 1; Rho = 2.62*10^-12 Exp[-(500000 - r[t] )/64500]; Rtet = Rearth + 500000 - r[t] ; ugc = 6.67*10^11; massEarth = 5.9721986*10^24; v = Sqrt[(ugc*massEarth)/Rtet] c1 = 2.67; a1 = 8; m1 = 1000; c2 = 2.5; a2 = 4; m2 = 500; B1 = c1*a1/m1; B2 = c2*a2/m2; hh = 64500; dt = 0.005; k1 = 1; k2 = 1; (*Disturbance*) Q1 = v^2*r[t]*Rho*Cos[Theta]Cos[Phi]*((Mu*(B1-B2)/2)- (dt*hh)*(1-Cos^2[Phi[t]]Sin^2[Theta[t]])^0.5/Cos[Theta[t]]Cos[Phi[t]]); Q2 = v^2*r[t]*Rho*Sin[Theta]Sin[Phi]*((Mu*(B2-B1)/2)- (dt*hh)*(1-Cos^2[Phi[t]]Sin^2[Theta[t]])^0.5/Cos[Theta[t]]Cos[Phi[t]]); (*Flip Flop Control Law*) u1 = k1*Theta'[t]-(Q2+u2)*(Theta'[t]+1)/Phi'[t]; u2 = k2*Phi'[t]-(Q1+u1)*Phi'[t]/(Theta'[t]+1); (*Anti Disturbance Control Law*) sigma = (u1+Q1/(2*Theta'[t]+1))+(u2+Q2/(2*Phi'[t])); (*System Dynamics*) e1 = Theta''[t]+(Theta'[t]+1)*(2*sigma- 2*Phi'[t]*Tan[Phi[t]])+3*Sin[Theta[t]]*Cos[Theta[t]]-Q1; e2 = Phi''[t]+2*sigma*Phi'[t]+ ((Theta'[t]+1)^2+3*Cos^2[Theta[t]])*Sin[Phi[t]]*Cos[Phi[t]]-Q2; e3 = r'[t]-sigma*r[t]; (*Initial conditions*) system = NDSolve[{e1 == 0, e2 == 0, e3 == 0, Theta == 0.1, Theta' == 0, Phi == 0.1, Phi' == 0.001, r == 1}, {Theta,Phi,r}, {t,0,tend}, MaxSteps -> Infinity]; (*Plot*) Plot[Evaluate[Theta[t] /.system], {t,0,tend}, Frame -> True, ImageSize -> {400, 300}, PlotStyle -> Red, LabelStyle -> Directive, FrameLabel -> {"time(s)", "Theta (rad)"}, PlotRange -> All] Plot[Evaluate[Phi[t] /.system], {t,0,tend}, Frame -> True, ImageSize -> {400, 300}, PlotStyle -> Blue, LabelStyle -> Directive, FrameLabel -> {"time(s)", "Phi (rad)"}, PlotRange -> All]  So the problem here lies in the (Flip Flop Control Law ) section. Since the equation in u1 is dependent on u2 and vice-versa, therefore this results in an algebraic looping and one method I know of in order to overcome this problem is to introduce a unit delay. In SIMULINK, the "Unit Delay" block is available, but in Mathematica, how do I program a unit delay for this system so I can simulate and plot the results successfully for Theta and Phi ? Thank you so much for your effort and time. I truly appreciate it and attached is a Mathematica program file to this problem for your reference. Best regards, Aaron Aw Attachments: