Unable to perform assignment because the left and right sides have a different number of elements.

% differential equation (SWING EQUATION)
%d^2(y)/d^2(t) = (-a*D)/(2*H)*dy/dt + (a*P)/(2*H) -(G*B*a)/(2*h*z)*sin(y-c);
%%MATLAB CODE%%
clc;
clear all;
N=100000;
T=linspace(0,100,N);
H=T(2)-T(1);
c=0.1; % theta B
H1=2.37; %Inertia
D= 0.002: 0.016; %damping term
P=0.85; %power
a=1; %constant angular velocity
G=1; %voltage of bus
f=60; %excitation frequesncy
s1=1; %psi v
B=1 + 0.1*cos(f*T + s1); %voltage of bus
z=0.645; %transient reactance
F=@(X2)X2;
G=@(X1,X2,T)-(a.*D)./(2.*H1).*X2 + (a.*P)./(2.*H1) -(G.*B.*a)./(2.*H1.*z).*((X1 - (((X1).^3)./6)).*cos(c) -(1-(((X1).^2)./2).*sin(c)));
X1=zeros(1,N);
X2=zeros(1,N);
X1(1)=2;
X2(1)=0;
for i=1:N-1
    K1=F(X2(i));
    L1=G(X1(i),X2(i),T(i));
    K2=F(X2(i)+((H/2)*L1));
    L2=G(X1(i)+((H/2)*K1),X2(i)+((H/2)*L1),T(i));
    K3=F(X2(i)+((H/2)*L2));
    L3=G(X1(i)+((H/2)*K2),X2(i)+((H/2)*L2),T(i));
    K4=F(X2(i)+((H)*L3));
    L4=G(X1(i)+((H)*K3),X2(i)+((H)*L3),T(i));
    X1(i+1)=X1(i)+(H/6)*(K1+2*K2+2*K3+K4);
    X2(i+1)=X2(i)+(H/6)*(L1+2*L2+2*L3+L4);
end
figure(1)
plot(T,X1)
hold on
plot(T,X2,'r');
legend('x_1(t)','x_2(t)');
xlabel('Normalized Time');ylabel('States of the system')
figure(2)
plot(X1,X2)
xlabel('x_1');ylabel('x_2')

Results:

Unable to perform assignment because the left and right sides have a different number of elements.

Error in TRIAL (line 31)

X1(i+1)=X1(i)+(H/6)*(K1+2*K2+2*K3+K4);

Can someone please help me with this code? I do not know where I went wrong. It is the swing eqaution and i need to sove using rungekutta method. The output i get is also methioned above!

[KSSV]

Your L1, L2, L3, L4 all are 1X100000 vectors..so this line:

X1(i+1)=X1(i)+(H/6)*(K1+2*K2+2*K3+K4);

RHS gives you a 1X100000 vector which you cannot save as single element.

[Bhairavi]

How do I corrrect it please help! Thank you

[KSSV]

B should be an array or a scalar?

B=1 + 0.1*cos(f*T + s1); %voltage of bus

[Bhairavi]

scalar

[Bhairavi]

I need to get figures like this and I am stuck

[VBBV]

B=1 + 0.1*cos(f*T + s1); %voltage of bus

In your definition it is vector with size of 1x100000. Try with a scale as below , it should work.

t = 1;% e.g. 
B=1 + 0.1*cos(f*t + s1); %voltage of bus

[KSSV]

I suspect B should be having a time difference?

B=1 + 0.1*cos(f*dT + s1); %voltage of bus

If you take T, it is a vector, and where veer B is used, it will be a vector. Check the formula for B.

Answer 1

[James Tursa]

Why are you replacing the sin(y-c) term with approximations? Why haven't you just programmed it as is? It appears you have done this:

sin(y-c) = sin(y)*cos(c) - cos(y)*sin(c)

and then for some reason replaced the sin(y) and cos(y) with truncated Taylor Series approximations, e.g.

sin(y) ~ y - y^3/6

cos(y) ~ 1 - y^2/2

Well, if you know that y is going to remain very small then maybe you can get away with this, but I frankly don't see the point of doing it in the first place unless this eventually will be hosted on a system with no trig functions available.

And for B, since it depends on time, you need to calculate it for the specific time involved in a calculation and make it a part of the inputs to the function handles that depend on it inside the loop. E.g.,

G=@(X1,X2,T,B) etc.  % <-- Added B to the input argument list

Although your equations as written don't appear to explicitly depend on T, you are not adjusting T in your RK4 scheme and B depends on this.

Finally, you have two different G's in your code. One is a double variable and the other is a function handle. I suppose this will work since the G in the function handle body gets used prior to the function handle being assigned the same name. But it is confusing code at best. I would advise changing the name of one of those G's. They are both annotated with the same comment "voltage of bus". I suspect that "1" in the B formula is the G=1 from above? If so, use the (changed) variable name rather than a "magic" number. Similar comment for the 0.1 in the B formula.

An outline of what needs to happen with your code is:

G = @(X1,X2,T,B) -(a.*D)./(2.*H1).*X2 + (a.*P)./(2.*H1) -(G.*B.*a)./(2.*H1.*z).*((X1 - (((X1).^3)./6)).*cos(c) -(1-(((X1).^2)./2).*sin(c))); % added B input
B = @(T) 1 + 0.1*cos(f*T + s1); % voltage of bus a function of time
    :
for i=1:N-1
    K1=F(X2(i));
    L1=G(X1(i),X2(i),T(i),B(T(i))); % altered for time
    K2=F(X2(i)+((H/2)*L1));
    L2=G(X1(i)+((H/2)*K1),X2(i)+((H/2)*L1),T(i)+H/2,B(T(i)+H/2)); % altered for time
    K3=F(X2(i)+((H/2)*L2));
    L3=G(X1(i)+((H/2)*K2),X2(i)+((H/2)*L2),T(i)+H/2,B(T(i)+H/2)); % altered for time
    K4=F(X2(i)+((H)*L3));
    L4=G(X1(i)+((H)*K3),X2(i)+((H)*L3),T(i)+H,B(T(i)+H)); % altered for time
    X1(i+1)=X1(i)+(H/6)*(K1+2*K2+2*K3+K4);
    X2(i+1)=X2(i)+(H/6)*(L1+2*L2+2*L3+L4);
end

You could cut down on some of this code by using a state vector instead of individual variables.