Issues with ode45
2 ビュー (過去 30 日間)
古いコメントを表示
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
0 件のコメント
採用された回答
Star Strider
2021 年 9 月 29 日
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.
.
4 件のコメント
Star Strider
2021 年 9 月 29 日
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.
.
その他の回答 (0 件)
参考
製品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!