plotting this program is causing a problem how to fix it?

2 views (last 30 days)
clear
k = 1E-3;
[y(4),y(5),y(6)] = meshgrid(-4:0.25:4,-4:0.25:4,-4:0.25:4);
Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.
dy(1) = (P - y(1).*((abs(y(4)))^2 +1))./tf;
dy(2) = (P - y(2).*((abs(y(5)))^2 +1))./tf;
dy(3) = (P - y(3).*((abs(y(6)))^2 +1))./tf;
dy(4)= (y(1)-a).*((y(4))./tc) + (k./tc).*(y(5)).*cos(y(7));
dy(5)= (y(2)-a).*((y(5))./tc) + (k./tc).*(y(4)).*cos(y(7)) + (k./tc).*(y(6))*cos(y(8));
dy(6)= (y(3)-a).*((y(6))./tc) + (k./tc).*(y(5)).*cos(y(8));
dy(7) = o(1,1) - (k./tc).*((y(4)./y(5)) + (y(5)./y(4))).*sin(y(7)) + (k./tc).*(y(6)/y(5)).*sin(y(8));
dy(8) = o(1,2) - (k./tc).*((y(5)./y(6)) + (y(6)./y(5))).*sin(y(8)) + (k./tc).*(y(4)/y(5)).*sin(y(7));
r = sqrt((dy(4)).^2 + (dy(5)).^2 + (dy(6)).^2 );
quiver(y(4),y(5),y(6), (dy(4))./r, (dy(5))./r, (dy(6))./r,1/2,'linewidth',1);
hold on;
axisa equal
set(gca,'Fontsize',20)
axis([-4 4 -4 4 -4 4])
while (true)
y0 = ginput(1);
ti = 0;
tf = 1E-2;
tspan=[ti tf];
[T,Y]= ode45(@(t,y) rate_eq(t,y,k),tspan,y0);
plot3(Y(:,4),Y(:,5),Y(:,6));
xlabel("A1");
ylabel("A2");
zlabel("A3");
end
function dy = rate_eq(t,y,k)
dy = zeros(8,1);
P = 0.2;
a = 0.1;
tf = 230E-6;
tc = 30E-9;
o(1,1) = 3E4;
o(1,2) = 4E4;
dy(1) = (P - y(1).*((abs(y(4)))^2 +1))./tf;
dy(2) = (P - y(2).*((abs(y(5)))^2 +1))./tf;
dy(3) = (P - y(3).*((abs(y(6)))^2 +1))./tf;
dy(4)= (y(1)-a).*((y(4))./tc) + (k./tc).*(y(5)).*cos(y(7));
dy(5)= (y(2)-a).*((y(5))./tc) + (k./tc).*(y(4)).*cos(y(7)) + (k./tc).*(y(6))*cos(y(8));
dy(6)= (y(3)-a).*((y(6))./tc) + (k./tc).*(y(5)).*cos(y(8));
dy(7) = o(1,1) - (k./tc).*((y(4)./y(5)) + (y(5)./y(4))).*sin(y(7)) + (k./tc).*(y(6)/y(5)).*sin(y(8));
dy(8) = o(1,2) - (k./tc).*((y(5)./y(6)) + (y(6)./y(5))).*sin(y(8)) + (k./tc).*(y(4)/y(5)).*sin(y(7));
end
% is my algorithm is correct in this? please modify if thereis any possible scenario.

Accepted Answer

Walter Roberson
Walter Roberson on 5 Oct 2022
[y(4),y(5),y(6)] = meshgrid(-4:0.25:4,-4:0.25:4,-4:0.25:4);
You want y(4), y(5), y(6) to each represent a matrix of values. You cannot do that with () indexing, but you could do
[y{4},y{5},y{6}] = meshgrid(-4:0.25:4,-4:0.25:4,-4:0.25:4);
Then you have
dy(1) = (P - y(1).*((abs(y(4)))^2 +1))./tf;
ON the right you would have to change that y(4) to y{4}, like
(P - y(1).*((abs(y{4}))^2 +1))./tf;
and then you would have to adjust the ^2 for the fact that you are working with matrices
(P - y(1).*((abs(y{4})).^2 + 1))./tf;
but you have not defined y(1) at this point. If you did use [y{4},y{5},y{6}] = meshgrid(etc) above then y(1) would have been created implicity, but it would be an empty cell and so
dy(1) = (P - y{1}.*((abs(y{4})).^2 + 1))./tf;
would be calculating with empty y{1} and so the calculation is going to return an empty array as a result, so you would further need
dy{1} = (P - y{1}.*((abs(y{4})).^2 + 1))./tf;
but that is a complicated way to calculate what is going to be dy(1) = {[]};
Well, except for the fact that P is not defined at that point, and tf is not defined either.
dy(4)= (y(1)-a).*((y(4))./tc) + (k./tc).*(y(5)).*cos(y(7));
a and tc and k are not defined either... and you have not assigned anything to y(7) ...
dy(7) = o(1,1) - (k./tc).*((y(4)./y(5)) + (y(5)./y(4))).*sin(y(7)) + (k./tc).*(y(6)/y(5)).*sin(y(8));
o is not defined either...
It looks to me as if you are trying to create a bunch of formulas for differential equations. If you were to syms various variables into existence and use syms y [1 7] then you could create the formuals. But you would not be able to plot the formulas until you substituted in numeric values for the symbolic variables.
I recommend that you have a look at the Symbolic Toolbox, in particular at the work-flow used in the first example for odeFunction
  1 Comment
SAHIL SAHOO
SAHIL SAHOO on 5 Oct 2022
I tried and there was lot of chaos I will leave it here though and thanks for help.

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by