Divergent solution using ode45

I have a code (see below) that solve the second order differential equation. As a output it gives the position and velocity as a function of time. For some values of omega 0.48,0.5,0.85 it gives the diverging solution of position. I guess this is due to some error but I can not find the error. How to overcome this errors for above omegas. Please help regarding this. I am new to MATLAB.

%% MAIN PROGRAM %%
t=linspace(0,1000,5000);
y0=[1 0];
omega=0.85;
[tsol, ysol]=ode45(@(t,y0) firstodefun4(t,y0,omega), t, y0, omega);
position=ysol(:,1);
velocity=ysol(:,2);

%% FUNCTION DEFINITION%%
function dy=firstodefun4(t,y0,omega)
F=1;gamma=0.01;omega0=1;
kappa=0.15;
dy=zeros(2,1);
dy(1)=y0(2);
dy(2)=2*F*sin(omega*t)-2*gamma*y0(2)-omega0^2*(1+4*kappa*sin(2*omega*t))*y0(1);
end

Answer 1

[Torsten]

We don't know the reason - technically, your code is correct. Although I would change it as shown below.

What makes you think that the solution you get is incorrect ? Do you have any indication that it should be bounded as t -> 00 ? For me, the behaviour looks very regular.

%% MAIN PROGRAM %%
t=linspace(0,50,5000);
y0=[1 0];
omega=0.85;
F=1;
gamma=0.01;
omega0=1;
kappa=0.15;
[tsol, ysol]=ode45(@(t,y0) firstodefun4(t,y0,omega,F,gamma,omega0,kappa), t, y0);
position=ysol(:,1);
velocity=ysol(:,2);
plot(tsol,position)

%% FUNCTION DEFINITION%%
function dy=firstodefun4(t,y0,omega,F,gamma,omega0,kappa)
  dy=zeros(2,1);
  dy(1)=y0(2);
  dy(2)=2*F*sin(omega*t)-2*gamma*y0(2)-omega0^2*(1+4*kappa*sin(2*omega*t))*y0(1);
end

[Abhik Saha]

Thanks for the suggestion regarding the code. The behaviour should not be divergent because for other cases of omega it gives sustained oscillating behaviour with non divergent behaviour this is reason I think that for this case it gives wrong result. Also if I take the average over time oscillation then it increases if the time span is high because it increases with increase in time.

[Torsten]

If your code above were wrong, it were wrong for all values of omega.

So the reason for the behaviour cannot be your code, but must either be the type of equation or the integrator (or both).

Concerning the equation, it can be such that small inaccuracies in the numerical solution can lead to big effects. Something you cannot really influence.

Concerning the integrator, you can try to strengthen the tolerances RelTol and AbsTol in order to damp the described effect as best as possible. Or try another integrator (e.g. RADAU5 from Hairer-Wanner).

Did you try to get a symbolic solution with kappa = 0 ? This can at least give you an impression of the exactness of the numerical solution.

syms t y(t)
eqn = diff(y,2)-2*sin(0.85*t)+0.02*diff(y)+y==0;
Dy = diff(y);
conds = [y(0)==1,Dy(0)==0];
y = dsolve(eqn,conds)

y = 

ynum = matlabFunction(y);
t=0:1:200;
hold on
plot(t,ynum(t))

y0=[1 0];
omega=0.85;
F=1;
gamma=0.01;
omega0=1;
kappa=0.0;
[tsol, ysol]=ode45(@(t,y0) firstodefun4(t,y0,omega,F,gamma,omega0,kappa), t, y0);
position=ysol(:,1);
velocity=ysol(:,2);
plot(tsol,position)
hold off

%% FUNCTION DEFINITION%%
function dy=firstodefun4(t,y0,omega,F,gamma,omega0,kappa)
  dy=zeros(2,1);
  dy(1)=y0(2);
  dy(2)=2*F*sin(omega*t)-2*gamma*y0(2)-omega0^2*(1+4*kappa*sin(2*omega*t))*y0(1);
end

So ode45 performs well in this case.

[Abhik Saha]

I have checked with RelTol and Abstol but it does not change much I have also check with smaller time steps but again the result remains same. If possible can you try with RADAU5 from Hairer-Wanner because I do not know how to use this. HJust to check whether it gives divergent result or not. Apart from that I have not tried with kappa=0. But I have tried with other ode integrator it does not change much.

[Jano]

Is it possible for this to accept a range of value for omega, for example omega = [0.85:5].

Answer 2

[Sam Chak]

Hi @Abhik Saha

If you reduce the value for , then the system should be stable.

For more info, please check out Floquet theory.

% MAIN PROGRAM %
n      = 400;
t      = linspace(0, n*100, n*10000+1);
y0     = [1 0];
omega  = 0.85;
F      = 1;
gamma  = 0.01;
omega0 = 1;
kappa  = 0.1469;

[tsol, ysol] = ode45(@(t, y) firstodefun4(t, y, omega, F, gamma, omega0, kappa), t, y0);

position = ysol(:, 1);
velocity = ysol(:, 2);
plot(tsol, position)

% ODE
function dy = firstodefun4(t, y, omega, F, gamma, omega0, kappa)
    dy    = zeros(2,1);
    dy(1) = y(2);
    dy(2) = 2*F*sin(omega*t) - 2*gamma*y(2) - omega0^2*(1 + 4*kappa*sin(2*omega*t))*y(1);
end

[Abhik Saha]

Hi @Sam Chak

Thank you for the comment. But I have checked that still there are few range of omegas like 0.47-0.50, 0.86-1.1 where I get the divergent solution. For those cases reducing the kappa value does not change the result it shows the divergent behaviour. Could you please check how to deal with above mentioned omegas to get sustained oscillation?

[Sam Chak]

@Abhik Saha, I changed the value of κ (via trial-and-error) because I thought of reducing the effect of the periodic function in the term. The original should be restored. You can definitely run the simulations for different values of angular frequency ω. There is nothing wrong about the response being divergent as (unstable).

You need to work out the algebra in order to find out fundamental matrix solution of the system, so that you can determine the stability codition following from the Floquet theorem. You can watch some videos here:

https://www.youtube.com/watch?v=N_zmYDnjACs

https://www.youtube.com/watch?v=3TbVWzRY2hU

by Prof. Nathan Kutz.

[Abhik Saha]

@Sam Chak Thank you for sharing the information related Floquet theorem. I will watch the videos and if I find some criteria then I will post here.

[Nikoo]

@Sam Chak

I was implementing the unstability and divergence boundries of a system with 2dof in MATLAB, but my matrices have periodic coefficient.

clc
clear all ;

% Define parameters 

N = 2 ; %Number of blades

I_thetaoverI_b = 2 ; % Moment of inertia pitch axis over I_b 
I_psioverI_b = 2 ; % Moment of inertia yaw axis over I_b 

C_thetaoverI_b = 0.00; % Damping coefficient over I_b 
C_psioverI_b = 0.00; % Damping coefficient over I_b 

h = 0.3; % rotor mast height, wing tip spar to rotor hub
hoverR = 0.34;
R = h / hoverR;

gamma = 4; % lock number

V = 325 ; % the rotor forward velocity [knots]
Omega = V/R; % the rotor rotational speed [RPM]

freq_1_over_Omega = 1 / Omega;

%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;

%the propeller aerodynamic coefficients
H_u = V/2;

%%%%%%%%%%%the flap moment aerodynamic coefficients for small V
%M_b = -1*(1 + V^2)/8 ;
%M_u = V/4;
%the propeller aerodynamic coefficients
%H_u = (V^2/2)*log(2/V);

f_pitch= 0.01:5:140;
f_yaw= 0.01:5:140;

phi_steps = linspace(0, pi, 50); % Evaluation points from 0 to pi

divergence_map = [];

Rdivergence_map = [];

unstable = [];

for i = 1:length(f_pitch)
    for j = 1:length(f_yaw)
        for phi = phi_steps

            % Calculate stiffness for the current frequency
            w_omega_pitch = 2*pi*f_pitch(i);
            w_omega_yaw = 2*pi*f_yaw(j);
        
            K_psi = (w_omega_pitch^2) * I_psioverI_b;
            K_theta = (w_omega_yaw^2) * I_thetaoverI_b;
        
        
            % Define inertia matrix [M]
            M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
                         -sin(2*phi), I_psioverI_b + 1 - cos(2*phi)];
            
            % Define damping matrix [D]
            D11 =  h^2*gamma*H_u*(1 - cos(2*phi)) - gamma*M_b*(1 + cos(2*phi)) - (2 + 2*h*gamma*M_u)*sin(2*phi);
            D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) - 2*(1 + cos(2*phi)) - 2*h*gamma*M_u*cos(2*phi);
            D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 - cos(2*phi)) - 2*h*gamma*M_u*cos(2*phi);
            D22 =  h^2*gamma*H_u*(1 + cos(2*phi)) - gamma*M_b*(1 - cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);
            
            D_matrix = [D11, D12;
                        D21, D22];
            
            % Define stiffness matrix [K]
            K11 = K_theta - h*gamma*V*H_u*(1 - cos(2*phi)) + gamma*V*M_u*sin(2*phi);
            K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
            K21 = -h*gamma*V*H_u*sin(2*phi) - gamma*V*M_u*(1 - cos(2*phi));
            K22 = K_psi - h*gamma*V*H_u*(1 + cos(2*phi)) - gamma*V*M_u*sin(2*phi);
            
            K_matrix = [K11, K12;
                        K21, K22];
        
            A_top = [zeros(2, 2), eye(2)];
            A_bottom = [-inv(M_matrix) * K_matrix, -inv(M_matrix) * D_matrix];
            A = [A_top; A_bottom];
        
            eigenvalues = eig(A);        
        
            % Stability condition 
            % Flutter
            
            if any(real(eigenvalues) > 0)
                unstable = [unstable; K_psi, K_theta];
            end
            
            % Divergence condition
            if det(K_matrix) < 0
                divergence_map = [divergence_map; K_psi, K_theta];
            end 
        
            % 1/Ω *(Divergence) proximity check
            for ev = eigenvalues'
                if abs(ev - freq_1_over_Omega) < 1e-2
                    Rdivergence_map = [Rdivergence_map; K_psi, K_theta];
                end
            end
        end 
    end
end

% Plot the Flutter and divergence maps
figure;
hold on;

scatter(unstable(:,1), unstable(:,2), 'filled');
scatter(divergence_map(:,1), divergence_map(:,2), 'filled', 'r');
scatter(Rdivergence_map(:, 1), Rdivergence_map(:, 2), 'filled', 'g');

xlabel('K\_psi');
ylabel('K\_theta');

title('Whirl Flutter Diagram');
legend('Flutter area', 'Divergence area', ' 1/Ω Divergence area');

hold off;

Could you take a look at my approach? The final plot isn't implemented correctly.

[Sam Chak]

Hi @Nikoo

Some engineering stuff here as I can see the eigenvalues of a two coupled 2nd-order system. I would suggest you to post as New Question to attract more attentions from the Aeroelastic flutter experts. If possible, attach a sketch of the expected result or plot.