Swinging spring ODE system in Matlab - How do I make a position vector follow the path?

Question

I'm very new to Matlab. I have a script using ode45 and arrow.m to show the motion of a swinging spring with a mass on Matlab as it moves through 3-D space. The program is almost doing what I want it to. Right now the density of the diamonds is showing the speed of the spring in effect (except for when ode45 takes it personal favorite number samples) and the velocity is almost accurately accounted for with the step size of the function (at least with the speed that my computer is running the code). What I want to do with this is have the position vector that I commented out in the code only show up at the instantaneous position of the mass, that is the endpoint of the curve, and not at every point that the diamonds show up at. I've looked around for help, but everything I try seems to cause errors. If somebody could point me in the right direction it would be much appreciated. Please try running the program and you will see what I mean, also play around with the parameters of the function.

function elasticPendulum(totalTime)
time=1;
totalTime=20
X=1
Vx=0
Y=0
Vy=0
Z=1.1
Vz=0
w=3;
g=9;
l=1;
while (time<totalTime)
    tspan=[0,time];
    x0=[X,Vx,Y,Vy,Z,Vz];
    [t,x]=ode45(@F,tspan,x0);
    plot3(x(:,1),x(:,3),x(:,5),'dk'), axis([-2 2 -2 2 -10 2]);
    grid on;
    o=[0, 0, 0];
    Xeq=[0, 0, -(g/(w^2)+l)];
    arrow(o,Xeq,'Length',5);
    %{ 
    Xt=[x(:,1) x(:,3) x(:,5)]
    arrow(o,Xt,'Length',5);
    %}
    time=time+.1*(((x(2))^2+(x(4))^2+(x(6))^2)^(1/2))
    pause(.1);
end

    function xprime=F(t,x)
    r=sqrt(x(1)^2+x(3)^2+x(5)^2);
    xprime=[x(2);-w*(r-l)/r*x(1);x(4);-w*(r-l)/r*x(3);x(6);-w*(r-l)/r*x(5)-g];
    end

end

Answer 1

I think that you can accomplish what you want simply by setting Xt to this so that only the last vector is plotted on each iteration:

Xt = [x(end,1) x(end,3) x(end,5)];

Also, you mention that things are working "except for when ode45 takes it personal favorite number samples." You can can tell ode45 to use fixed step output simply by changing tspan:

tspan = 0:dt:time;

where dt = 1/time (or you could use linspace). This is not the same thing as using a fixed step solver – read my answer to this question to maybe gain some understanding on this.

I don't think that you're updating your animation properly. You update time, but not the initial conditions. Therefore, you integrate over the same path on each iteration of the while loop. Is there a reason that you're adjusting the end time on each iteration? Are you trying to find the oscillation period or something?

Also, your method of animation is rather crude/inefficient. You should read about the output options that can be set via odeset for Matlab's ODE solvers. In particular 'OutputFcn'. You might even be able to use and/or look at the code for some of the built in output functions – for example, type edit odeplot or edit odephas3 in your command window. Here is a simple output function I created for your case:

function elasticPendulum
totalTime = 20;
stepsPerSec = 10;
X = 1; Vx = 0;
Y = 0; Vy = 0;
Z = 1.1; Vz = 0;
w = 3; g = 9; l = 1;

opts = odeset('OutputFcn',@(t,x,flag)arrowplot(t,x,flag,w,l,g));
tspan = linspace(0,totalTime,totalTime*stepsPerSec);
x0 = [X,Vx,Y,Vy,Z,Vz];
ode45(@(t,x)F(t,x,w,l,g),tspan,x0,opts);

function xprime=F(t,x,w,l,g)    %#ok<INUSL>
r=sqrt(x(1)^2+x(3)^2+x(5)^2);
xprime=[x(2);-w*(r-l)/r*x(1);x(4);-w*(r-l)/r*x(3);x(6);-w*(r-l)/r*x(5)-g];

function status=arrowplot(t,x,flag,w,l,g)   %#ok<INUSL>
persistent h; % Handle for moving arrow
status = 0;
o = [0, 0, 0];
switch flag
    case 'init'
        axis([-2 2 -2 2 -10 2]); grid on; hold on;
        Xeq = [0, 0, -(g/(w^2)+l)];
        arrow(o,Xeq,'Length',5);
        plot3(x(1,:),x(3,:),x(5,:),'dk');  % Initial position
        Xt = [x(1,end) x(3,end) x(5,end)];
        h = arrow(o,Xt,'Length',5);        % Initial arrow, get handle
    case 'done'
        hold off; status = 1;
    otherwise
        plot3(x(1,:),x(3,:),x(5,:),'dk');        % Plot new positions
        Xt = [x(1,end) x(3,end) x(5,end)];
        arrow(h,'Start',o,'Stop',Xt,'Length',5); % Update arrow
        pause(0.1);
end

The call to plot3 on every iteration (the otherwise case in arrowplot) is still inefficient as it adds new high level plot objects, taking up more memory. The better/faster way is to get a handle from the first call to plot3 and use set and get to add the new data points. See the code for odeplot and odephas3 for how you might do this (it's a bit advanced).

Note how I also pass parameters via anonymous functions rather than by creating sub-functions. This is a bit of a matter of style.