I'm not getting any output here, can anyone please tell me what's the issue here?
Afficher commentaires plus anciens
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)
Réponses (1)
Torsten
le 30 Déc 2022
for i=1:(length(time)-1)
instead of
for i=1:(length(t)-1)
Catégories
En savoir plus sur Ordinary Differential Equations dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
