Continuous Errors Trying to Solve for an 11-State Variable w/ ODE45

2 views (last 30 days)
I am trying to solve for 11 different variables by defining a function that includes 11 ODEs, 1 for each variable. These ODEs are coupled so they must be solved simultaneously. I have gone over my code a bunch of times and cannot find any errors, but when I run it I keep getting an error that says my functino returns a vector of length 8 which is not compatible with my 11 initial conditions. Please let me know if anyone knows how to fix this! Here is my code:
% Parameter Values (rad/s)
a1=0.51;
a2=0.251;
a3=0.151;
omega1=0.182;
omega2=0.424;
omega3=0.827;
%X=[qBarc1,qBarc2,qBarc3,qBarc4,qBar1,qBar2,qBar3,qBar4,w1,w2,w3]';
X0=[0.2000;0.1000;0.4472;0.86600;0.1035;0.0518;0.2315;0.9659;0;0;0];
odes=@(t,X)[(1/2)*(-a3*sin(omega3*t)*X(2)+a2*sin(omega2*t)*X(3))+(1/2)*X(4)*a1*sin(omega1*t); % 1st ODE,qBarc1
(1/2)*(a3*sin(omega3*t)*X(1)-a1*sin(omega1*t)*X(3))+(1/2)*X(4)*a2*sin(omega2*t); % 2nd ODE,qBarc2
(1/2)*(-a2*sin(omega2*t)*X(1)+a1*sin(omega1*t)*X(2))+(1/2)*X(4)*a3*sin(omega3*t); % 3rd ODE,qBarc3
-(1/2)*(a1*sin(omega1*t)*X(1)+a2*sin(omega2*t)*X(2)+a3*sin(omega3*t)*X(3)); % 4th ODE,qBarc4
(1/2)*(-X(11)*X(6)+X(10)*X(7))+(1/2)*X(8)*X(9); % 5th ODE,qBar1
(1/2)*(X(11)*X(5)-X(9)*X(7))+(1/2)*X(8)*X(10); % 6th ODE,qBar2
(1/2)*(-X(10)*X(5)+X(9)*X(6))+(1/2)*X(8)*X(11); % 7th ODE,qBar3
-(1/2)*(X(9)*X(5)+X(10)*X(6)+X(11)*X(7)); % 8th ODE,qBar4
-(-1200*diff(a1*sin(omega1*t)) + 100*(diff(X(10))-diff(a2*sin(omega2*t))) - 200*(diff(X(11))-diff(a3*sin(omega3*t))) ...
+ 639.68*(X(9)-a1*sin(omega1*t)) + 44.422*(X(10)-a2*sin(omega2*t)) - 71.075*(X(11)-a3*sin(omega3*t)) ...
+ 341.09*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8)) ...
+ 19.739*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))...
-25.266*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/1200; % 9th ODE,w1
-(100*(diff(X(9))-diff(a1*sin(omega1*t))) - 2200*diff(a2*sin(omega2*t)) + 300*(diff(X(11))-diff(a3*sin(omega3*t))) ...
+ 53.307*(X(9)-a1*sin(omega1*t)) + 977.29*(X(10)-a2*sin(omega2*t)) + 106.61*(X(11)-a3*sin(omega3*t)) ...
+ 28.424*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8)) ...
+ 434.26*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))...
+37.899*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/2200; % 10th ODE,w2
-(-200*(diff(X(9))-diff(a1*sin(omega1*t))) + 300*(diff(X(10))-diff(a2*sin(omega2*t))) - 3100*diff(a3*sin(omega3*t)) ...
- 106.61*(X(9)-a1*sin(omega1*t)) + 133.27*(X(10)-a2*sin(omega2*t)) + 1101.7*(X(11)-a3*sin(omega3*t)) ...
- 56.849*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8)) ...
+ 59.218*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))...
+ 391.63*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/3100]; % 11th ODE,w3
[t,XSol]=ode45(odes,[0 60],X0);
Error using odearguments
@(T,X)[(1/2)*(-A3*SIN(OMEGA3*T)*X(2)+A2*SIN(OMEGA2*T)*X(3))+(1/2)*X(4)*A1*SIN(OMEGA1*T);(1/2)*(A3*SIN(OMEGA3*T)*X(1)-A1*SIN(OMEGA1*T)*X(3))+(1/2)*X(4)*A2*SIN(OMEGA2*T);(1/2)*(-A2*SIN(OMEGA2*T)*X(1)+A1*SIN(OMEGA1*T)*X(2))+(1/2)*X(4)*A3*SIN(OMEGA3*T);-(1/2)*(A1*SIN(OMEGA1*T)*X(1)+A2*SIN(OMEGA2*T)*X(2)+A3*SIN(OMEGA3*T)*X(3));(1/2)*(-X(11)*X(6)+X(10)*X(7))+(1/2)*X(8)*X(9);(1/2)*(X(11)*X(5)-X(9)*X(7))+(1/2)*X(8)*X(10);(1/2)*(-X(10)*X(5)+X(9)*X(6))+(1/2)*X(8)*X(11);-(1/2)*(X(9)*X(5)+X(10)*X(6)+X(11)*X(7));-(-1200*DIFF(A1*SIN(OMEGA1*T))+100*(DIFF(X(10))-DIFF(A2*SIN(OMEGA2*T)))-200*(DIFF(X(11))-DIFF(A3*SIN(OMEGA3*T)))+639.68*(X(9)-A1*SIN(OMEGA1*T))+44.422*(X(10)-A2*SIN(OMEGA2*T))-71.075*(X(11)-A3*SIN(OMEGA3*T))+341.09*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8))+19.739*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))-25.266*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/1200;-(100*(DIFF(X(9))-DIFF(A1*SIN(OMEGA1*T)))-2200*DIFF(A2*SIN(OMEGA2*T))+300*(DIFF(X(11))-DIFF(A3*SIN(OMEGA3*T)))+53.307*(X(9)-A1*SIN(OMEGA1*T))+977.29*(X(10)-A2*SIN(OMEGA2*T))+106.61*(X(11)-A3*SIN(OMEGA3*T))+28.424*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8))+434.26*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))+37.899*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/2200;-(-200*(DIFF(X(9))-DIFF(A1*SIN(OMEGA1*T)))+300*(DIFF(X(10))-DIFF(A2*SIN(OMEGA2*T)))-3100*DIFF(A3*SIN(OMEGA3*T))-106.61*(X(9)-A1*SIN(OMEGA1*T))+133.27*(X(10)-A2*SIN(OMEGA2*T))+1101.7*(X(11)-A3*SIN(OMEGA3*T))-56.849*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8))+59.218*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))+391.63*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/3100]
returns a vector of length 8, but the length of initial conditions vector is 11. The vector returned by
@(T,X)[(1/2)*(-A3*SIN(OMEGA3*T)*X(2)+A2*SIN(OMEGA2*T)*X(3))+(1/2)*X(4)*A1*SIN(OMEGA1*T);(1/2)*(A3*SIN(OMEGA3*T)*X(1)-A1*SIN(OMEGA1*T)*X(3))+(1/2)*X(4)*A2*SIN(OMEGA2*T);(1/2)*(-A2*SIN(OMEGA2*T)*X(1)+A1*SIN(OMEGA1*T)*X(2))+(1/2)*X(4)*A3*SIN(OMEGA3*T);-(1/2)*(A1*SIN(OMEGA1*T)*X(1)+A2*SIN(OMEGA2*T)*X(2)+A3*SIN(OMEGA3*T)*X(3));(1/2)*(-X(11)*X(6)+X(10)*X(7))+(1/2)*X(8)*X(9);(1/2)*(X(11)*X(5)-X(9)*X(7))+(1/2)*X(8)*X(10);(1/2)*(-X(10)*X(5)+X(9)*X(6))+(1/2)*X(8)*X(11);-(1/2)*(X(9)*X(5)+X(10)*X(6)+X(11)*X(7));-(-1200*DIFF(A1*SIN(OMEGA1*T))+100*(DIFF(X(10))-DIFF(A2*SIN(OMEGA2*T)))-200*(DIFF(X(11))-DIFF(A3*SIN(OMEGA3*T)))+639.68*(X(9)-A1*SIN(OMEGA1*T))+44.422*(X(10)-A2*SIN(OMEGA2*T))-71.075*(X(11)-A3*SIN(OMEGA3*T))+341.09*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8))+19.739*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))-25.266*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/1200;-(100*(DIFF(X(9))-DIFF(A1*SIN(OMEGA1*T)))-2200*DIFF(A2*SIN(OMEGA2*T))+300*(DIFF(X(11))-DIFF(A3*SIN(OMEGA3*T)))+53.307*(X(9)-A1*SIN(OMEGA1*T))+977.29*(X(10)-A2*SIN(OMEGA2*T))+106.61*(X(11)-A3*SIN(OMEGA3*T))+28.424*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8))+434.26*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))+37.899*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/2200;-(-200*(DIFF(X(9))-DIFF(A1*SIN(OMEGA1*T)))+300*(DIFF(X(10))-DIFF(A2*SIN(OMEGA2*T)))-3100*DIFF(A3*SIN(OMEGA3*T))-106.61*(X(9)-A1*SIN(OMEGA1*T))+133.27*(X(10)-A2*SIN(OMEGA2*T))+1101.7*(X(11)-A3*SIN(OMEGA3*T))-56.849*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8))+59.218*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))+391.63*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/3100]
and the initial conditions vector must have the same number of elements.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
I have double checked the equations and they should all be correct. Each equation is equal to the derivative of the variable associated with that state. Please let me know if anything is wrong or why this only produces 8 variables when there are 11 equations.

Answers (1)

Walter Roberson
Walter Roberson on 6 May 2022
-(-1200*diff(a1*sin(omega1*t))
omega1 is a scalar constant. At any instant in an ode45() call, t is a scalar value. omega1*t is therefore a scalar value, and sin() of it would be a scalar value. Multiply by the scalar a1 to get a scalar value.
diff() of a scalar value is empty.
Reminder:
diff() has two very different uses. When applied to a symbolic function or symbolic expression, diff() does a calculus differentiation. When applied to a numeric scalar or array, diff() is numeric differences, similar to diff(X) --> X(:,2:end) - X(:,1:end-1)
If you want the derivative of a1*sin(omega1*t) with respect to time, and you are not working symbolically, then do the devivative manually at the time you write the formula, getting a1*cos(omega1*t)
  2 Comments
Erin Hayes
Erin Hayes on 6 May 2022
Thank you so much, I did not realize that that would return a scalar. I will try this!
Erin Hayes
Erin Hayes on 6 May 2022
I just replaced all the diff() that had t in them with their derivatives. This is my new code:
% Parameter Values (rad/s)
a1=0.51;
a2=0.251;
a3=0.151;
omega1=0.182;
omega2=0.424;
omega3=0.827;
%X=[qBarc1,qBarc2,qBarc3,qBarc4,qBar1,qBar2,qBar3,qBar4,w1,w2,w3]';
X0=[0.2000;0.1000;0.4472;0.86600;0.1035;0.0518;0.2315;0.9659;0;0;0];
odes=@(t,X)[(1/2)*(-a3*sin(omega3*t)*X(2)+a2*sin(omega2*t)*X(3))+(1/2)*X(4)*a1*sin(omega1*t); % 1st ODE,qBarc1
(1/2)*(a3*sin(omega3*t)*X(1)-a1*sin(omega1*t)*X(3))+(1/2)*X(4)*a2*sin(omega2*t); % 2nd ODE,qBarc2
(1/2)*(-a2*sin(omega2*t)*X(1)+a1*sin(omega1*t)*X(2))+(1/2)*X(4)*a3*sin(omega3*t); % 3rd ODE,qBarc3
-(1/2)*(a1*sin(omega1*t)*X(1)+a2*sin(omega2*t)*X(2)+a3*sin(omega3*t)*X(3)); % 4th ODE,qBarc4
(1/2)*(-X(11)*X(6)+X(10)*X(7))+(1/2)*X(8)*X(9); % 5th ODE,qBar1
(1/2)*(X(11)*X(5)-X(9)*X(7))+(1/2)*X(8)*X(10); % 6th ODE,qBar2
(1/2)*(-X(10)*X(5)+X(9)*X(6))+(1/2)*X(8)*X(11); % 7th ODE,qBar3
-(1/2)*(X(9)*X(5)+X(10)*X(6)+X(11)*X(7)); % 8th ODE,qBar4
-(-1200*a1*omega1*cos(omega1*t) + 100*(diff(X(10))-a2*omega2*cos(omega2*t)) - 200*(diff(X(11))-a3*omega3*cos(omega3*t)) ...
+ 639.68*(X(9)-a1*sin(omega1*t)) + 44.422*(X(10)-a2*sin(omega2*t)) - 71.075*(X(11)-a3*sin(omega3*t)) ...
+ 341.09*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8)) ...
+ 19.739*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))...
-25.266*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/1200; % 9th ODE,w1
-(100*(diff(X(9))-a1*omega1*cos(omega1*t)) - 2200*a2*omega2*cos(omega2*t) + 300*(diff(X(11))-a3*omega3*cos(omega3*t)) ...
+ 53.307*(X(9)-a1*sin(omega1*t)) + 977.29*(X(10)-a2*sin(omega2*t)) + 106.61*(X(11)-a3*sin(omega3*t)) ...
+ 28.424*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8)) ...
+ 434.26*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))...
+37.899*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/2200; % 10th ODE,w2
-(-200*(diff(X(9))-a1*omega1*cos(omega1*t)) + 300*(diff(X(10))-a2*omega2*cos(omega2*t)) - 3100*a3*omega3*cos(omega3*t) ...
- 106.61*(X(9)-a1*sin(omega1*t)) + 133.27*(X(10)-a2*sin(omega2*t)) + 1101.7*(X(11)-a3*sin(omega3*t)) ...
- 56.849*(X(4)*X(5)+X(3)*X(6)-X(2)*X(7)-X(1)*X(8)) ...
+ 59.218*(-X(3)*X(5)+X(4)*X(6)+X(1)*X(7)-X(2)*X(8))...
+ 391.63*(X(2)*X(5)-X(1)*X(6)+X(4)*X(7)-X(3)*X(8)))/3100]; % 11th ODE,w3
[t,XSol]=ode45(odes,[0 60],X0);
I got the same error as before even with replacing these. Could the error be coming from when I use diff() on the state variables such as diff(X(9))? Is there another way to indicate derivative? Thanks again for your help!

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!