While loops in ode45
7 views (last 30 days)
Show older comments
José Anchieta de Jesus Filho
on 11 May 2021
Commented: José Anchieta de Jesus Filho
on 17 May 2021
Hello guys, I would like to ask for your help. Well, I'm trying to find the amplitude function using ode45, but I'm not having success at the end of the code. I am submitting the exact solution and the solution that I am implementing in matlab.
Basically, my goal is to be able to extract the amplitudes of the steady-state solution from ode45 and plot according to each frequency. To do this, I started with the loop that will do the calculation for each frequency and, within that loop, it must verify that the average of the peaks minus the last peak is less than the error, if this error is greater, it must continue computing the loop.
When this error is less, I would like to take all of these peaks and capture the last peak, so "amp (i) = y1 (end);" and, at the end of it, he computes each amp as a function of w. I would like help to solve my problem.
Only one detail, apparently this code works for the first frequency, but not for others. Then, he can find exactly the amplitude that the exact equation, at the initial frequency.
I0 = 0.086;
k1 = 4.9085e+03;
c1 = 0.001*k1;
l = 0.4665;
a = 0.1841;
F1 = 1.3688;% N
w = 2:0.5:14; %Hz
% IC
theta0 = 0; omega0 = 0;
IC2 = [theta0,omega0];
% Time
t0 = 0;
tf = 1;
maxerr = 1e-4;
erro = inf;
%exact solution
d1 = (k1*a^2 - I0*(w*2*pi).^2).^2 + (c1*(a^2)*w*2*pi).^2;
X = F1*l./sqrt(d1);
for i = 1:length(w)
while erro >= maxerr
tf = tf + 1;
tt = [t0 tf];
f1 = @(time) F1.*l.*sin(w(i)*time*2*pi);
%sdot
sdot2 = @(t,x) [x(2);
(f1(t) - (c1.*a^2).*x(2) - (k1*a^2)*x(1))/I0];
% numerical integration
[time,state_values2] = ode45(sdot2,tt,IC2);
theta1 = state_values2(:,1);
% maxim
y1 = findpeaks(theta1);
erro = abs(y1(end) - mean(y1));
end
amp(i) = y1(end);
end
disp(X(1))
disp(erro)
disp(tf)
disp(amp)
figure(2)
plot(w,X,'--k',w,amp,'ok');
l1 = ' Matlab Solution';
l2 = ' Analytical solution';
legend(l1,l2);
axis square
xlim([w(1) 14])
ylim([0 0.1])
ax1 = gca;
ax1.YAxis.Exponent = -2;
set(gca,'FontName','Times New Roman')
set(gca,'FontSize',20)
set(gca,'linewidth',1.5)
xlabel('\omega (Hz)','FontName','Cambria Math','FontSize',20)
ylabel('X (rad)','FontName','Cambria Math','FontSize',20)
set(gcf,'color','w');
box off
0 Comments
Accepted Answer
Jan
on 11 May 2021
Edited: Jan
on 11 May 2021
You do not reset erro to inf after the first iteration. Then all following iterations are not entered.
Replace:
erro = inf;
...
for i = 1:length(w)
while erro >= maxerr
by
...
for i = 1:length(w)
erro = inf;
while erro >= maxerr
Your code wastes a lot of time with repeated calculations. Instead of starting the integration at t0 repeatedly, you can start at the former tf and use the previous final value as new initial value.
for i = 1:length(w)
f1 = @(time) F1.*l.*sin(w(i)*time*2*pi);
sdot2 = @(t,x) [x(2);
(f1(t) - (c1.*a^2).*x(2) - (k1*a^2)*x(1))/I0];
t0 = 0;
tf = 1;
IC = IC2;
y1 = [];
erro = inf;
while erro >= maxerr
tt = [t0, t0 + 1]; % 1 additional second
[time,state_values2] = ode45(sdot2, tt, IC2);
theta1 = state_values2(:, 1);
y1 = [y1, findpeaks(theta1)]; % Append new peaks
erro = abs(y1(end) - mean(y1));
% New initial values:
t0 = tt(2);
IC2 = state_values2(end, :);
end
amp(i) = y1(end);
end
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!