The initial subscript problem was the result of:
y0=0;
because there are 12 differential equations. There are several typographical errors that I did my best to correct, replacing ‘f’ with ‘fr’, ‘q’ with ‘y’, and ‘A’ with ‘A1’. Even with those changes and expanding ‘y0’, it will not run in the 55 seconds allowed here.
I encourage you to experiment with it —
P.ph1_l1=-0.012247303909990;
P.ph1_l2=0.635309628411612;
P.rm=1.509423275703383;
P.rs=1.960632785814804;
P.Cp=1.170700000000000e-07;
P.ohm=33.720923694042130;
P.A1=0.089447240667428;
P.p1=9.769092596233520;
P.Ca=0.001400000000000;
P.rf=0.5;
P.B=0.028248587570621;
P.Io=0.554857887911763;
P.Ia=-0.002864751824000;
P.Ib=0.952971893800000;
P.pc1=1.140600000000000;
P.Cs=0.010100000000000;
P.niu=9.520415341114659e-07;
P.fm=2;
P.phi_lm=0.503431090400000;
P.I11=0.2;
P.Zeeta=3.696152028525914e-04;
P.C1=0.839958085300000;
P.rM=4.241322704123775;
P.N=-0.006529965080390;
P.G1=1.363350182000000;
P.Rl=0.0002;
tspan=0:0.001:100;
% y0=0;
y0 = zeros(12,1);
opts = odeset('Mass',@(t,q) mass(t,q,P));
[t,q] = ode15s(@(t,q) fu(t,q,P),tspan,y0,opts);
figure
plot(t,q)
grid
legend(compose('y_%d', 1:12), 'Location','bestoutside')
function H = mass(t,y,P)
ph1_l1=P.ph1_l1;
ph1_l2=P.ph1_l2;
rm=P.rm;
rs=P.rs;
Cp=P.Cp;
ohm=P.ohm;
A1=P.A1;
p1=P.p1;
Ca=P.Ca;
rf=P.rf;
B=P.B;
Io=P.Io;
Ia=P.Ia;
Ib=P.Ib;
pc1=P.pc1;
Cs=P.Cs;
niu=P.niu;
fm=P.fm;
phi_lm=P.phi_lm;
I11=P.I11;
Zeeta=P.Zeeta;
C1=P.C1;
rM=P.rM;
N=P.N;
G1=P.G1;
Rl=P.Rl;
H = zeros(12,12);
H(1,1)=1;
H(2,2)=1;
H(2,8)=ph1_l1*rm;
H(2,10)=ph1_l2*rm;
H(3,3)=1;
H(4,4)=A1*rm/rs;
H(4,8)=rm/rs;
H(5,5)=1;
H(6,6)=A1*rm/rs;
H(6,10)=rm/rs;
H(7,7)=1;
H(8,8)=1;
H(9,9)=1;
H(10,10)=1;
H(11,11)=Cp*ohm;
H(12,12)=Cp*ohm;
end
function dydt = fu(t,y,P)
ph1_l1=P.ph1_l1;
ph1_l2=P.ph1_l2;
rm=P.rm;
rs=P.rs;
Cp=P.Cp;
ohm=P.ohm;
A=P.A1;
p1=P.p1;
Ca=P.Ca;
rf=P.rf;
B=P.B;
Io=P.Io;
Ia=P.Ia;
Ib=P.Ib;
pc1=P.pc1;
Cs=P.Cs;
niu=P.niu;
fm=P.fm;
phi_lm=P.phi_lm;
I11=P.I11;
Zeeta=P.Zeeta;
C1=P.C1;
rM=P.rM;
N=P.N;
G1=P.G1;
Rl=P.Rl;
dydt = [y(2)
-p1^4*y(1)-Ca*y(2)+ph1_l1*rm*rf^2*B*sin(rf*t)+ph1_l2*rm*rf.^2*B*sin(rf*t)+(rf.^2*B*sin(rf*t-rf*cos(rf*t))*Io)
y(4)
-A*y(3)*pc1^4-(Cs*pc1^4*A+Ca*A)*y(4)-Ca*y(8)*I11-y(11)*niu*(Ia-Ib)+fm*phi_lm+((rm/rs)*rf^2*B*sin(rf*t)-Ca*rf*B*cos(rf*t))*I11
y(6)
-A*y(5)*pc1^4-(Cs*pc1^4*A+Ca*A)*y(6)-Ca*y(10)*I11-y(12)*niu*(Ia-Ib)+fm*phi_lm+((rm/rs)*rf^2*B*sin(rf*t)-Ca*rf*B*cos(rf*t))*I11
y(8)
rf^2*ph1_l1*y(1)+Zeeta*ph1_l1*y(2)+2*C1*(rs/rM)*pc1^3*y(3)-rf^2*y(7)-Zeeta*y(8)+rf^2*B*sin(rf*t)
y(10)
rf^2*ph1_l2*y(1)+Zeeta*ph1_l2*y(2)+2*C1*(rs/rM)*pc1^3*y(5)-rf^2*y(9)-Zeeta*y(10)+rf^2*B*sin(rf*t)
-N*G1*y(4)-y(11)/Rl
-N*G1*y(6)-y(12)/Rl];
end
EDIT — (20 Nov 2023 at 14:20)
Changed integrator from ode45 to ode15s.
.
