How do I output the second derivative from an ODE solver for further use?

I am currently working on a coupled ode problem, and using ode45 solver to solve this.

My code is something like this:

bubble.m

function rdot = f(t, r)

%Some mathematical expressions
....
...
...

rdot(1) = r(2); % first derivative of r1
rdot(2) = (X11 - X21 - (X31*(X41 - X51)))/X81; % second derivative of r1

rdot(3) = r(4); % first derivative of r2
rdot(4) = (X12 - X22 - (X32*(X42 - X52)))/X82; % second derivative of r2

rdot = rdot';

And,

bubble_plotter.m

clc;
clear all;
%close all;
time_range = [0 1d-4];
r1_eq = 10d-6;
r2_eq = 1d-6;
f_s = 20000;
T_s = 1/f_s;
initial_conditions = [r1_eq 0.d0 r2_eq 0.d0];
[t, r] = ode45('bubble', time_range, initial_conditions);
r1_us=1000000*r(:,1);
r1dot=r(:,2);
r2_us=1000000*r(:,3);
r2dot=r(:,4);
normalized_time = t./T_s;
figure(101);
h1 = plot(normalized_time, r1_us, 'b-', normalized_time, r2_us, 'k-.');
set(h1,'linewidth',3);
txt = text(0.03, 21.5, (['d = ' num2str(d_close), '\mum']));
txt.FontSize = 25;
legend('Bubble1', 'Bubble2')
legend('Location','Northwest')
set(gca,'xscale','linear','FontSize',26)
set(gca,'yscale','linear','FontSize',26)
set(gca,'XMinorTick','on')
set(gca,'YMinorTick','on')

You can see, I am pulling out r(:,1), r(:,2), r(:,3) and r(:,4) as r1_us, r1dot and r2_us, r2dot for plotting and further uses. But I also need the values of rdot(2) and rdot(4) from the main function file. How can I pull those double derivative values of r1 and r2 into my plotting file for further use?

Answer 1

[Torsten]

After the line

[t, r] = ode45('bubble', time_range, initial_conditions);

insert

for i=1:numel(t)
  t_actual = t(i);
  r_actual = r(i,:);
  rdot(i,:) = bubble(t_actual,r_actual);
end

Best wishes

Torsten.

[Vikash Pandey]

@ Torsten, I am little confused. You can see that I can already retrieve first derivatives of r1 and r2 as:

r1dot=r(:,2);
r2dot=r(:,4);

I am interested in the double derivatives that are calculated as

rdot(2) = (X11 - X21 - (X31*(X41 - X51)))/X81; % second derivative of r1

rdot(4) = (X12 - X22 - (X32*(X42 - X52)))/X82; % second derivative of r2

The way you have written the expression

rdot(i,:) = bubble(t_actual,r_actual);

means that

rdot is calculated for i = 1 to length(t), for ":" that could correspond to 1 and 2. Is my understanding correct?

So if I wanted to use double derivative of r1 in my calculation, the variable name should be rdot(:,1) and for double derivative of r2 it should be rdot(:,2)? Correct?

[Vikash Pandey]

By the way actual r1 and r2 are r(:,1) and r(:,3) respectively.

[Torsten]

In my code, rdot(:,2) and rdot(:,4) give you the second derivatives of r(:,1) and r(:,3).

[Vikash Pandey]

I will check and let you know the results soon. By the way I also found an alternative to obtain differentiation of a curve through this:

r1dd0t = diff(r1dot)./diff(t);

Hope this is correct too.

[Torsten]

Yes, but less exact.

Why not use the explicit formulae for the second derivatives if you just supplied them in "bubble" ?

[Vikash Pandey]

@ Torsten, I inserted this following your advise.

for i=1:numel(t)
  t_actual = t(i);
  r_actual = r(i,:);
  rdot(i,:) = bubble_mettin(t_actual,r_actual);
end

r1ddot = rdot(:,2);
r2ddot = rdot(:,4);

And, I got good results, I believe. The result (thick blue bubble1 curve) matches with the implict (thin) curve that I obtained using "diff".

Hope I am finally correct. can you confirm, see the relevant section of the code carefully please for the second derivative.

[Torsten]

Yes, that's what I meant.