Issues with ode45

I'm simulating the dynamics for an autonomous system and and the system is described as follows: x1dot=-ax1+x2, x2dot=-bx2+u, udot=-cu. My initial conditions for [x1 x2 u]=[0 0 0]. I am randomly generaring values of x1, x2, and u for 100 iterations and also randomly generating values for the constants a, b, and c. My code is pasted below. I keep getting the error "Not enough input arguments". I also would like to have x1(t), x2(t), and u(t) on three separate plots not one. If anyone has any suggestions it is 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] = ode45(@(x1,x2,u) dynamics(x1,x2,u,a,b,c), tspan, x0, options );

plot(t,y), hold on

end

title('Problem 1 System Dynamics')

xlabel('Time')

ylabel('Dynamic System')

hold off

function system = dynamics(x1,x2,u,a,b,c)

system=[-a*x1+x2;-b*x2+u;-c*u] ;

end

Answer 1

[Star Strider]

There are two problems (and one suggestion to improve its efficiency):

First, ‘dymanics’ should only present the independent and dependent variables to ode45.

Second, ‘x2’ is a vector, and needs to be subscripted as such. Change the subscripts if my guess as to their assignments is not correct.

Third, it may be ‘stiff’, so I substituted ode15s for ode45, since the integration seemed to be taking longer than I expected it to.

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(@(x1,x2) dynamics(x1,x2,u,a,b,c), tspan, x0, options );
    plot(t,y), hold on
end
title('Problem 1 System Dynamics')
xlabel('Time')
ylabel('Dynamic System')
hold off

function system = dynamics(x1,x2,u,a,b,c)
system = zeros(3,1);
system(1) = -a*x1+x2(1);
system(2) = -b*x2(2)+u;
system(3) = -c*u;
% system=[-a*x1+x2;-b*x2+u;-c*u] 
end

Experiment to get different results.

.

[Matthew Tortorella]

This seems along the lines of what I expected, but why do you only subscript x2 and not x1. And why in the function is x2 subscripted as x2(1) in system(1) but x2 is subscripted as x2(2) in system(2). I think I may just not be understanding why x2 is a vector and x1 is not. I was referencing the following from a matlab syntax explanation and maybe this will help you understand my confusion:

y′1=y2

y2'=μ(1−y1^2)y2−y1.

function dydt = vdp1(t,y)

%VDP1 Evaluate the van der Pol ODEs for mu = 1

%

% See also ODE113, ODE23, ODE45.

% Jacek Kierzenka and Lawrence F. Shampine

% Copyright 1984-2014 The MathWorks, Inc.

dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

[Star Strider]

‘... but why do you only subscript x2 and not x1’

Since ‘x1’ is a scalar in each call to ‘dynamics’ there is not need to subscript it.

‘I think I may just not be understanding why x2 is a vector and x1 is not.’

Here, ‘x1’ is equivalent to a time variable, so only one scalar instance of it exists in each call to ‘dynamics’.

In the code you quoted, ‘y’ is the variable to be integrated. The time variable ‘t’ exists in the ‘vdp’ function, however it is not used. (It could be, however it just isn’t in that system.)

You probably want ‘dynamics’ instead to be:

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

With that change (and concordant changes in the ode15s call), the code produces this result —

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 );
    plot(t,y), hold on
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

It appears similar to the earlier result, however now the code is correct.

.

[Matthew Tortorella]

@Star Strider that does make more sense for x1 and x2. Thank you for clearing that up. Do you have any idea on how to plot the separate equations making up the dynamics (system(1), system(2), sysytem(3)) on their own plot when running the code a single time? This would allow me to see how x1(t), x2(t), and u(t) behave over the 100 iterations but how the code is now I am seeing all three of them on one plot.

[Star Strider]

My pleasure.

I was off doing other things most of today. As for the plots, I would use fewer iterations in the ‘ii’ loop (perhaps 10, rather than 100), and then use subplot plots to view the individual runs. Creating subplot plots in a loop is straightforward.

.