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

4 次查看(过去 30 天)
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?

采纳的回答

Torsten
Torsten 2018-11-14
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.
  7 个评论
Vikash Pandey
Vikash Pandey 2018-11-14
编辑:Vikash Pandey 2018-11-14
@ 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.

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Interactive Control and Callbacks 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by