ODE45 giving complex solution unexpectedly

Hi, I'm trying to solve the following odes, which are two coupled 2nd order odes giving a simple representation of a bearing. The equations are:

The functions I have are as shown:

function eqns_of_motion()
    
    close all
    clear
    clc
    
    % Define parameters
    M = 0.6; % Mass
    C = 200; % Damping coefficient
    %K = bearingStiffness(7.94,46.98,31.1,4.2082,4.2082,200,0.3); % Stiffness
    K = 8.3687e+09;
    % cr = 1; % Constant
    Wr = 100; % External force
    Z = 6; % Number of balls 
    Ns = 1500; % shaft speed in rpm
    Db = 7.94; % ball diameter
    Dc = 39.04; % pitch/cage diameter
    
    ws = 2*pi*Ns/60; % shaft anglular velocity
    wc = (ws/2)*(1-Db/Dc); % cage angular velocity

    %%%function definition, see odefunc below too for the summation term
    fx = @(t, x, y, dxdt, dydt) (-(C*dxdt)) / M;
    fy = @(t, x, y, dxdt, dydt) (-(C*dydt) - Wr) / M;

    start_time = 0;
    end_time = 0.7;
    no_steps = (end_time-start_time) / 10e-5;

        % Define initial conditions
    tspan = linspace(start_time,end_time,no_steps); % Time span
    initial_conditions = [10e-6; 0; 10e-6; 0]; % [x0; dxdt0; y0; dydt0]
    
   
    [t, u] = ode45(@(t, u) odefunc(t, u, fx, fy, M, K, Z, wc, ws), [0 0.7], initial_conditions);
    
    x = u(:, 1);
    dxdt = u(:, 2);
    y = u(:, 3);
    dydt = u(:, 4);
    
    % acceleration
    d2xdt2 = gradient(dxdt, t); 
    d2ydt2 = gradient(dydt, t); 
    
    figure;
    
    subplot(2,2,1);
    plot(t, x, 'r');
    xlabel('Time');
    ylabel('x(t)');
    title('Solution of x(t)');
    
    subplot(2,2,2);
    plot(t, y, 'b');
    xlabel('Time');
    ylabel('y(t)');
    title('Solution of y(t)');

    subplot(2,2,3);
    plot(t, d2xdt2, 'r');
    xlabel('Time');
    ylabel('d^2x/dt^2');
    title('Acceleration of x(t)');

    subplot(2,2,4);
    plot(t, d2ydt2, 'b');
    xlabel('Time');
    ylabel('d^2y/dt^2');
    title('Acceleration of y(t)');

    figure; 
    Fs = 1/(t(2)-t(1));
    L = length(x);
    f = Fs*(0:(L/2))/L;
    X = fft(x);
    Y = fft(y);

    P2x = abs(X/L);
    P2y = abs(Y/L);

    subplot(2,1,1);
    plot(f, P2x(1:L/2+1), 'r');
    title('frequency spectrum of x(t)');
    xlabel('Frequency (Hz)');
    ylabel('|P1(f)|');

    subplot(2,1,2);
    plot(f, P2y(1:L/2+1), 'b');
    title('frequency spectrum of y(t)');
    xlabel('Frequency (Hz)');
    ylabel('|P1(f)|');

end

function dudt = odefunc(t, u, fx, fy, M, K, Z, wc, ws)
    x = u(1);
    dxdt = u(2);
    y = u(3);
    dydt = u(4);

        % Compute the summation term
    sum_term_x = 0;
    sum_term_y = 0;
    for i = 1:Z
        thetai = wc*t+2*pi*(i-1)/Z;
        

        sum_term_x = sum_term_x + K * (x * cos(thetai) + y * cos(thetai))^1.5 * cos(thetai);
        sum_term_y = sum_term_y + K * (x * sin(thetai) + y * sin(thetai))^1.5 * sin(thetai);
        
    end
    
    dudt = zeros(4, 1);
    dudt(1) = dxdt;
    dudt(2) = fx(t, x, y, dxdt, dydt) - sum_term_x/M - ws^2*cos(ws*t);  
    dudt(3) = dydt;
    dudt(4) = fy(t, x, y, dxdt, dydt) - sum_term_y/M - ws^2*sin(ws*t); 
end

   

But the output of both x and y both have complex components as when I try to plot displacement and acceleration I get a warning stating that imaginary parts are being ignored.

I have tried to interegate the ode45 function which shows the imaginary parts are appearing quite early in the solution, which would suggest an error in setting up the odes to be solved (I think!).

Any advice on how to debug this would be greatly appreciated as my experience in MATLAB is quite limited, thank you

Answer 1

[John D'Errico]

This seems clear. The only way in your code you can generate complex numbers is by raising a negative number to a fractional power. That happens in two places.

(x * cos(thetai) + y * cos(thetai))^1.5

(x * sin(thetai) + y * sin(thetai))^1.5

So, if (x * cos(thetai) + y * cos(thetai)) or its close cousin ever go even slightly negative, things will start to go to hell. And sadly, complex numbers are like wire coat hangars in your closet. They multiply. (Hmm. That seems like sort of a pun, because complex numbers in fact do propagate when you multiply them. Not intentional though.)

My point is, look carefully at what is happening. Is the solution doing somethign meaningful, BEFORE the complex numbers arise? If it is, then consider if ODE45 is just hitting the mathematical equivalent of a brick wall at that point. It does not understand that beyond this point lie dragons.

Anyway, it MIGHT also be just that you have a problem in your code. Or, perhaps the approximation inherant in that sum is merely inadequate.

[Matteo Thornton]

Ahh thanks :) in all the time spent looking for issues in the actual code I'd completely forgot that ^1.5 involves taking the square root...

In which case I would imagine the issue lies in that. I'll have a deeper look into it.

Also to answer your question, I wouldn't say the solution is meaningful as it definitely isn't correct, for example the x and y displacements shouldn't be trending up/down respectively, and the accelerations seems to have

an overly high frequency as shown below:

But while it isn't correct, it looks like it is on the right track, so hopefully looking into the fractional power will shed some light on it.