differentiation of the ODE solution

Hallo guys, I had ODE in my code. I solved it using an ode15i function. My code is such that the inputs of the ode15i function change when the output of the ODE reaches a limit. I did successfully solve the ODE and got the plot. Now I want to differentiate the solution to get the accelerations. I tried using diff() operator and passed the numerical outputs of the ODE but the graph seems incorrect (seems discontinous with sudden peaks and jumps).

Is there any way that I can differentiate the numerical solution / plot obtained from the ODE .

Answer 1

[Star Strider]

There are several possibilities, depending on what your ODE function is. The easiest way is to give your integrated solutions to your ODE function, since it codes for all the derivatives.

For example (using an anonymous function for the ODE function):

ode_fcn = @(t,y) ...;                                       % ODE Function
[t,y] = ode15i(ode_fcn, ...);                               % Integrate ODE
dy = ode_fcn(t,y);                                          % Calculate Derivatives

That is how I calculate the derivatives.

[Torsten]

But t is vector and y is a matrix ...

ode_fcn usually expects t as a scalar and y as a vector.

Best wishes

Torsten.

[Star Strider]

Good point, Torsten, thank you. It’s been a while since I actually did that.

It is necessary to use a loop for each value of the ‘t’ vector and each row of the ‘y’ matrix.

The corrected code would be:

ode_fcn = @(t,y) ...;                                       % ODE Function
[t,y] = ode15i(ode_fcn, ...);                               % Integrate ODE
for k1 = 1:size(t,1)
    dy(k1,:) = ode_fcn(t(k1),y(k1,:));                      % Calculate Derivatives
end

The ‘t’ vector will be the same for the integrated ODE and the calculated derivatives.

NOTE — Since I do not have your code to experiment with, this is UNTESTED CODE.

[yashwant kolluru]

Strider! I already did that, the problem was it was outputting only the last value of the function into my main.

Anyway, I declared dy as global variable and stored it in my ode function. But the problem is that the plot of the derivatives is discontinuous. The solution is 2nd order continuous and hence, I assume the derivatives also to be continuous (at least piece wise continuous). I get the graph with jumps at points where my ode breaks and starts with new set of values.

[yashwant kolluru]

Thanks Torsten! for the inputs!

[Star Strider]

I would not use global variables. There are many reasons to not use them, the most important being that it is very difficult to debug code that uses them. You can also inadvertently change them in other parts of your code.

I don’t have you code to experiment with. The for loop as I coded it should give you the matrix of your derivatives. Remember to index ‘dy’ as I did, because if you do not, only the last value will be returned.

[yashwant kolluru]

Ok I will try and explain you in a clear way.

I built a dynamics system of the form M*ydd(t) +P*yd(t) + Q*y(t) = h,

where ydd, yd are second and first derivatives of y. M, P, Q are 3X3 matrices and x is 3X1 column .

In order to solve it, I remodified the equation as below.

function dydt = expl(t,y,Mdash,Qdash,Pdash,hdash,Sdash) dydt = [y(4:6);(M^(-1)*h)-(M^(-1)*Q*y(1:3))];

The first three elements of y corresponds to displacements and next threee to velocties. I successfully obtain complete sol for y

The first three elements of dydt are same as last three of y, i.e velocities and next threee to accelertaions.

Now my main function as ode definition as follows

[t_expl,y_expl,te_expl,ye_expl,ie_expl] = ode45(@(tspan,y)expl(tspan,w_init,Mdash,Qdash,Pdash,hdash,Sdash(o+1)),tspan,w_init,opts);

and the ode would break accoriding to the condtion given in the event function associated withe ode. And the new loop number with new values of M, K, P, h starts and the ode continues.

The below is the event function of the ode

function [value,isterminal,direction] = events_expl(t,y,Sdash)

% value = abs(y(2) - y(1))-Sdash;

if (abs(y(2) - y(1)) > Sdash)

value =0;

else

value =1;

end

isterminal =1;

direction =0;

Now the plots of y seem continuous and genuine.

But the plots of dydt as discontiuous and something false. In order to verify that dydt graphs were wrong.

I did plot the first three element columns of dydt. Theoretically it should be equal to plot of the last three element colums of y (see the definition of dydt in the ode function).

The graphs of dydt seems to have correct values(same as y) when runng on aloop. But when the loop number changes, the values of dydt have some jumps.

Hope I did not confuse you.

Thanks!

[yashwant kolluru]

I attached the 2 figs one generated by y(4:6) and 2nd from dydt(1:3). For better understanding.

[yashwant kolluru]

I solved it ! Apologize for the inconvineince anyways, thanks for the help

[Star Strider]

I don’t understand you differential equation, but if it’s producing the correct output, that may not be the issue.

It’s better to post .png images than .fig files. I’ve not looked at them.

If using your original ODE function to calculate the derivatives is not working (and it might not if you’re using an event function to create a discontinuous solution), The next option is to use the gradient function on each column of the ‘y’ matrix separately. However this requires that you produce regularly-spaced vectors to it, so your ‘tspan’ argument would have to be a vector, for example:

tspan = linspace(t_min, t_max, 50);

to create a 50-element regularly-spaced vector. Then use the gradient function to take the derivative of each column of ‘y’ separately. Do not give it ‘y’ as a matrix, because it will calculate the derivatives between the columns as well as along the columns.