How do I plot acceleration from my ODE45 function?

Hey guys, I'm doing a simple spring mass damper problem using ODE45 and I'd like to generate plots of the acceleration, along with the displacement and velocity plots I already have. This is my main code:

   [t,y]=ode45(@(t,y) func(t,y),[0,100],[10,0]);

    figure(1)
    plot(t,y(:,1));
    xlabel('Time (s)')
    ylabel('Displacment (m)')

    figure(2)
    plot(t,y(:,2));
    xlabel('Time (s)')
    ylabel('velocity (m/s)')

   %figure(3)
   %plot(t,ydot(:,2))
   %xlabel('Time (s)')
   %ylabel('Acceleration (m/s/s)')

With this function being called:

function ydot=func(t,y)

    k=1000;
    c=1200;
    m=7000;

    ydot(1)=y(2);
    ydot(2)=((-k*y(1))/m) + ((-c*y(2))/m);
    ydot=ydot';

    end

What I'd like to plot is ydot(2) versus time but it doesn't stay in the work space when I run it without using the debugger and I get an error that tells me my ydot is not defined.

Any help would be appreciated.

Thanks, Conor.

Answer 1

[Jan]

Modify the function to be integrated such, that it accepts a matrix as input y:

function ydot=func(t, y)
k = 1000;
c = 1200;
m = 7000;
ydot(1, :) = y(2, :);
ydot(2, :) = -k * y(1, :) / m - c * y(2, :) / m;
end

By the way: Too many parentheses do not improve the clarity.

Then you can call func() with the output of the integration:

[t, y] = ode45(@func, [0,100], [10,0]);
ydot  = func(t, y);

Now you have ydot to the same times and positions as the integration steps.

By the way: @(t,y) func(t,y) creates an anonymous function, which calls func() with its original arguments. Then a normal function handle @func is cheaper and easier.

 but it doesn't stay in the work space

Correct. Each function has its own workspace and it is the purpose of functions, that they keep their internally used variables for themselves.

[Conor Downes]

Thanks for your help Jan!

Can you explain to me why ydot reruns only a 2x2 matrix when it's a function of t and y which are much larger matrices? This means i can't plot ydot versus t.

[Jan]

@Conor: Do you know the debugger already? See MATLAB: Debugger. It allows you to step through your code line by line. This would help you to find the mistake I made: ode45 provides the trajectory y as [numel(t) x 2] matrix, but I assumed a [2 x numel(t)] array. Simply transpose the matrix:

[t,y] = ode45(@(t,y) func(t,y),[0,100],[10,0]);
ydot  = func(t, y.').';  % <== Transpose twice
plot(t, y, t, ydot)

Now func() gets a [2 x numel(t)] array as input and for plotting the output ydot is transposed again to equal the [numel(t) x 2] orientation of y.

[Conor Downes]

Ah yes I see that now. Thanks.

[Daniel Gonzalez]

Can someone explain why this mathematically works? I'm having trouble seeing why the result of the second use of ode45 with the results from the first use is the acceleration.

[Jan]

Which 2nd use? func replies the derivative. ODE45 use it for an integration, but it you want the values of the derivative, you can use the output of func directly also. Here the outputs t and y of ODE45 are used to get the same time points and positions as during the integration.

There is no deeper mathematical meaning behing this operation.

[Daniel Gonzalez]

My apologies, I commented in the wrong thread.