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

3 次查看（过去 30 天）
Vikash Pandey
Vikash Pandey 2018-11-13
评论： Torsten 2018-11-14
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 个评论
显示 5更早的评论
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 Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by