While loops in ode45

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

Answer 1

[Jan]

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

[José Anchieta de Jesus Filho]

I would like to thank you for your response and for the time you took to answer me. I was able to complete my code successfully and I owe it to you. Thank you very much.

for i = 1:length(w)
    f1 = @(time) F1.*l.*sin(w(i)*time*2*pi);
    tf = 0;
    erro = inf;
    
    while erro >= maxerr          
            
        tf = tf + 10;    
        tt = [t0 tf];
        
        % 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,opts);    
        theta1 = state_values2(:,1);   
                
        % maxim   
        y1 = findpeaks(theta1);
        erro = abs((mean(y1) - y1(end))/y1(end));
                       
    end

    amp(i) = y1(end);
    
end