Trouble with 'global' & 'ode45' with plot in TDM spr-mass-damp. system

1 Ansicht (letzte 30 Tage)
Thomas Marshall
Thomas Marshall am 5 Dez. 2019
Bearbeitet: Walter Roberson am 5 Dez. 2019
I am having trouble trying to plot a response based on a Vibrations problem of a spring-mass-damper TDM system using the 'ode45' solver and 'global' commands (as asked by my instructor). My code below keeps giving me an error when trying to run the main .m file.
MAIN .m FILE
global ktmd
for i = 1:200
ktmd = 816000*i/100;
[t,x] = ode45(@ydot,(0:1),[0 0 0 0]);
z = x(((3/4)*end:end),1);
xmax(i,1) = max(z);
end
plot(ktmd*(1:200)/100,max)
FUNCTION FILE
function dy = ydot(t,y)
global ktmd
dy(1,1) = y(2);
dy(2,1) = -(816000/908)*y(1)-(ktmd/908)*y(1)-(32000/908)*y(2)+(ktmd/908)*y(3)+(125945.5*sin(29.98*t+0.86597));
dy(3,1) = y(4);
dy(4,1) = (ktmd/90.8)*y(1)-(ktmd/90.8)*y(3);
end
I get the error message about the "Index in position 1 is invalid" from the main .m file.
Does anyone know where i am going wrong with my code and trying to plot the response?
Thanks.

Antworten (1)

Walter Roberson
Walter Roberson am 5 Dez. 2019
Bearbeitet: Walter Roberson am 5 Dez. 2019
[t,x] = ode45(@ydot,(0:1),[0 0 0 0]);
z = x(((3/4)*end:end),1);
When you pass a vector of exactly 2 elements in as the timespan to ode45, then ode45 will return how-ever many rows of output that it feels are appropriate to achieve the integration tolerances. The number of rows it returns will often not be a multiple of 4, so 3/4*end will often not be an integer. You should use something like
z = x(ceil((3/4)*end):end,1);
Then on
plot(ktmd*(1:200)/100,max)
the first part is going to be length 200, and the second part is a call to the function max() without passing in any input arguments, which is going to give an error. You are also ignoring the xmax that you carefully calculated. Your second parameter should probably be xmax .
But then ktmd will be left as the last value it was assigned, which was from
ktmd = 816000*i/100;
with i = 200, and it would be that value that you would use as your plotting. It is much more likely that you would want your plot to be
plot(816000*(1:200)/100, xmax)

Kategorien

Mehr zu Plot Customization finden Sie in Help Center und File Exchange

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by