While loops in ode45
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
José Anchieta de Jesus Filho
el 11 de Mayo de 2021
Comentada: José Anchieta de Jesus Filho
el 17 de Mayo de 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 comentarios
Respuesta aceptada
Jan
el 11 de Mayo de 2021
Editada: Jan
el 11 de Mayo de 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
1 comentario
Más respuestas (0)
Ver también
Categorías
Más información sobre Annotations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!