# Issues plotting ode in a for loop

2 views (last 30 days)
Matthew Tortorella on 30 Sep 2021
I am simulating the dynamics of an autonomous system which is made up of following equations: x1dot=-ax1+x2, x2dot=-bx2+u, udot=-cu. The variables a,b,c are constants with applied constraints on them, the initial conditions for [x1 x2 u]=[0 0 0]. My code is posted below. I am trying to plot x1(t), x2(t), and u(t) on three separate plots for a given time span and for 100 random iterations of x1,x2,u,a,b,c. The issue I am running into is that my plots keep generating over and over for each iteration instead of placing the 100 trajectories for x1(t), x2(t), and u(t) on three separate plots in which all of the iterations for the given time span are on a single plot. Any help on this would be much appreciated. Thanks.
low=1;
high=1000;
ASa=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant a
ASb=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant b
ASc=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant c
ASx1=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable x1
ASx2=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable x2
ASu=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable u
for ii = 1:100 %signifies 100 iterations for each variable x1,x2,u,a,b,c
a =ASa(ii);
b=ASb(ii);
b = max([b, 1/a]); % Make sure b > 1/a
c=ASc(ii);
c = max([c, 1/(2*b)]); % Make sure c is > 1/(2b)
x1=ASx1(ii);
x2=ASx2(ii);
u=ASu(ii);
tspan = [0 5];
options = odeset('reltol', 1e-5, 'abstol', 1e-8 );
x0 = [0 0 0];
[t, y] = ode15s(@(t,x) dynamics(t,x,u,a,b,c), tspan, x0, options );
figure(1)
plot(t,y(:,1))
hold on
figure(2)
plot(t,y(:,2))
figure(3)
plot(t,y(:,3))
end
title('Problem 1 System Dynamics')
xlabel('Time')
ylabel('Dynamic System')
hold off
function system = dynamics(t,x,u,a,b,c)
system = zeros(3,1);
system(1) = -a*x(1)+x(2);
system(2) = -b*x(2)+u;
system(3) = -c*u;
% system=[-a*x1+x2;-b*x2+u;-c*u]
end

Sulaymon Eshkabilov on 30 Sep 2021
There are a couple of "hold on" s are missing and thus, not all simulated data get plotted. Here is the corrected part of the code:
low=1;
high=1000;
ASa=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant a
ASb=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant b
ASc=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for constant c
ASx1=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable x1
ASx2=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable x2
ASu=(high-low).*rand(100,1) + low;%column vector of uniformly distributed decimals between 1 and 1000 for variable u
for ii = 1:100 %signifies 100 iterations for each variable x1,x2,u,a,b,c
a =ASa(ii);
b=ASb(ii);
b = max([b, 1/a]); % Make sure b > 1/a
c=ASc(ii);
c = max([c, 1/(2*b)]); % Make sure c is > 1/(2b)
x1=ASx1(ii);
x2=ASx2(ii);
u=ASu(ii);
tspan = [0 5];
options = odeset('reltol', 1e-5, 'abstol', 1e-8 );
x0 = [0 0 0];
[t, y] = ode15s(@(t,x) dynamics(t,x,u,a,b,c), tspan, x0, options );
figure(1)
plot(t,y(:,1))
hold on
figure(2)
plot(t,y(:,2))
hold on % Necessary
figure(3)
plot(t,y(:,3))
hold on % Necessary
end
title('Problem 1 System Dynamics')
xlabel('Time')
ylabel('Dynamic System')
hold off
function system = dynamics(t,x,u,a,b,c)
system = zeros(3,1);
system(1) = -a*x(1)+x(2);
system(2) = -b*x(2)+u;
system(3) = -c*u;
% system=[-a*x1+x2;-b*x2+u;-c*u]
end
You may also consider to compute 1st and then plot the data, e.g.:
...
for ii = 1:100 %signifies 100 iterations for each variable x1,x2,u,a,b,c
....
[t, y] = ode15s(@(t,x) dynamics(t,x,u,a,b,c), tspan, x0, options );
T{ii}=t;
Y{ii}=y;
end
...
Then plot the results:
for ii=1:100
figure(1)
plot(T{ii},Y{ii}(:,1))
hold on
figure(2)
plot(T{ii},Y{ii}(:,2))
hold on
figure(3)
plot(T{ii},Y{ii}(:,3))
hold on
end
...