I'm not getting any output here, can anyone please tell me what's the issue here?
2 次查看(过去 30 天)
显示 更早的评论
h=0.01;
time=0:h:1;
a1(1)=rand(1,1);
a2(1)=rand(1,1);
g1(1)=rand(1,1);
g2(1)=rand(1,1);
o1(1)=((-3.14).*rand(1,1) + (3.14).*rand(1,1));
o2(1)=((-3.14).*rand(1,1) + (3.14).*rand(1,1));
s1(1)=((-3.14).*rand(1,1) + (3.14).*rand(1,1));
P = 0.5;
a = 0;
T = 2000;
k = 1E-3;
omega1 = 1E5;
omega2 = 2E5;
tp = 1E-12;
p = @(a1,g1) (P - g1 - (1 + 2.*g1).*(((a1)).^2))./T ;
q = @(a2,g2) (P - g2 - (1 + 2.*g2).*(((a2)).^2))./T ;
r = @(g1,a1,o2,o1) (g1).*(a1) + k.*(a2.*cos(o2-o1));
s = @(g2,a2,o2,o1) (g2).*(a2) + k.*(a1.*cos(o1-o2));
t = @(g1,a1,a2,o2,o1) (a.*g1 + omega1.*tp + k.*((a2./a1).*sin(o2-o1)));
u = @(g2,a1,a2,o2,o1) (a.*g2 + omega1.*tp + k.*((a1./a2).*sin(o1-o2)));
v = @(g1,a1,a2,o2,o1) (a.*(g1-g2) + k.*((a2./a1).*sin(o2-o1)) - k.*((a1./a2).*sin(o1-o2)) + omega1-omega2);
for i=1:(length(t)-1)
k1=p(a1(i),g1(i));
l1=q(a2(i),g2(i))
m1=r(g1(i),a1(i),o2(i),o1(i));
n1=s(g2(i),a2(i),o2(i),o1(i));
p1=t(g1(i),a1(i),a2(i),o2(i),o1(i));
q1=u(g2(i),a1(i),a2(i),o2(i),o1(i));
r1=v(g1(i),a1(i),a2(i),o2(i),o1(i));
k2=p((a1(i)+(h/2)*k1),(g1(i)+(h/2)*k1));
l2=q((a2(i)+(h/2)*l1),(g2(i)+(h/2)*l1));
m2=r((g1(i)+(h/2)*k1),(a1(i)+(h/2)*k1),(o2(i)+(h/2)*n1),(o1(i)+(h/2)*m1));
n2=s((g2(i)+(h/2)*l1),(a2(i)+(h/2)*l1),(o2(i)+(h/2)*n1),(o1(i)+(h/2)*m1));
p2=t((g1(i)+(h/2)*k1),(a1(i)+(h/2)*k1),(a2(i)+(h/2)*l1),(o2(i)+(h/2)*n1),(o1(i)+(h/2)*m1));
q2=u((g2(i)+(h/2)*l1),(a2(i)+(h/2)*l1),(a1(i)+(h/2)*k1),(o2(i)+(h/2)*n1),(o1(i)+(h/2)*m1));
r2=v((g1(i)+(h/2)*l1),(a1(i)+(h/2)*k1),(a2(i)+(h/2)*l1),(o2(i)+(h/2)*n1),(o1(i)+(h/2)*m1));
k3=p((a1(i)+(h/2)*k2),(g1(i)+(h/2)*k2));
l3=q((a2(i)+(h/2)*l2),(g2(i)+(h/2)*l2));
m3=r((g1(i)+(h/2)*k2),(a1(i)+(h/2)*k2),(o2(i)+(h/2)*n2),(o1(i)+(h/2)*m2));
n3=s((g2(i)+(h/2)*l2),(a2(i)+(h/2)*l2),(o2(i)+(h/2)*n2),(o1(i)+(h/2)*m2));
p3=t((g1(i)+(h/2)*k2),(a1(i)+(h/2)*k2),(a2(i)+(h/2)*l2),(o2(i)+(h/2)*n2),(o1(i)+(h/2)*m2));
q3=u((g2(i)+(h/2)*l2),(a2(i)+(h/2)*l2),(a1(i)+(h/2)*k2),(o2(i)+(h/2)*n2),(o1(i)+(h/2)*m2));
r3=v((g1(i)+(h/2)*l2),(a1(i)+(h/2)*k2),(a2(i)+(h/2)*l2),(o2(i)+(h/2)*n2),(o1(i)+(h/2)*m2));
k4=p((a1(i)+(h)*k3),(g1(i)+(h)*k3));
l4=q((a2(i)+(h)*l3),(g2(i)+(h)*l3));
m4=r((g1(i)+(h)*k3),(a1(i)+(h)*k3),(o2(i)+(h)*n3),(o1(i)+(h)*m3));
n4=s((g2(i)+(h)*l3),(a2(i)+(h)*l3),(o2(i)+(h)*n3),(o1(i)+(h)*m3));
p4=t((g1(i)+(h)*k3),(a1(i)+(h)*k3),(a2(i)+(h)*l3),(o2(i)+(h)*n3),(o1(i)+(h)*m3));
q4=u((g2(i)+(h)*l3),(a2(i)+(h)*l3),(a1(i)+(h)*k3),(o2(i)+(h)*n3),(o1(i)+(h)*m3));
r4=v((g1(i)+(h)*l3),(a1(i)+(h)*k3),(a2(i)+(h)*l3),(o2(i)+(h)*n3),(o1(i)+(h)*m3));
a1(i+1) = a1(i) + h*(k1+2*k2+2*k3+k4)/6;
a2(i+1) = g1(i) + h*(l1+2*l2+2*l3+l4)/6;
g1(i+1) = g1(i) + h*(m1+2*m2+2*m3+m4)/6;
g2(i+1) = a2(i) + h*(n1+2*n2+2*n3+n4)/6;
o1(i+1) = o1(i) + h*(p1+2*p2+2*p3+p4)/6;
o2(i+1) = o2(i) + h*(q1+2*q2+2*q3+q4)/6;
s1(i+1) = s1(i) + h*(r1+2*r2+2*r3+r4)/6;
end
plot(time,a2)
0 个评论
回答(1 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!