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

1 次查看（过去 30 天）
Bhairavi
Bhairavi 2022-6-23
编辑： James Tursa 2022-6-27
% 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!
  7 个评论
显示 5更早的评论
VBBV
VBBV 2022-6-23
编辑：VBBV 2022-6-23
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
KSSV 2022-6-23
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

James Tursa
James Tursa 2022-6-27
编辑：James Tursa 2022-6-27
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.

类别

在 Help Center 和 File Exchange 中查找有关 Runge Kutta Methods 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by