How do I output the second derivative from an ODE solver for further use?
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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?
0 commentaires
Réponse acceptée
Torsten
le 14 Nov 2018
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 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Interactive Control and Callbacks dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!