Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN

The following code uses `rk4` to simulate the dynamics defined via `fe` function. However, I encountered the warning `Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. `.

I felt that the function defined in `fe`, sometimes the Numerator equals zero, thus creates singularity. But I don't know how to avoid it.

Should I change the function formula (I prefer not to), or should I revise the code, or change the solver?

Thank you in advance for any suggestion.

clear;
global G
tic
n=40;
A = ones(n) - eye(n); 
G = graph(A);

num_samples =1;
p.G = A
dim_G = size(p.G);
p.K = 1;
p.N = dim_G(1); 

nIters = 2000;
tBegin = 0;
tEnd = 100;

% initial condition
Init = -pi + 2*pi.*rand(p.N,num_samples);
dim_thetaInit = size(Init);
Init = Init(:,1);

[t,sol] = rk4(@fe,tBegin,tEnd,Init,nIters,p);

function td = fe(t,theta,p)
[theta_i,theta_j] = meshgrid(theta);
adj_matrix = p.G;
a = 4;
b = 1;
td= sum(adj_matrix'.*( ( b*cos(theta_i-theta_j) * (2*a^2*cos(theta_i-theta_j)*sin(theta_i-theta_j)-2*b^2*cos(theta_i-theta_j)*sin(theta_i-theta_j)) ) / ( 2*(b^2*(cos(theta_i-theta_j))^2+a^2*(sin(theta_i-theta_j))^2)^(3/2) ) + (b*sin(theta_i-theta_j))/(  sqrt( b^2*cos(theta_i-theta_j)^2 + a^2*sin(theta_i-theta_j)^2 )   )  ),1)';
end

function [T,Y] = rk4(f,a,b,ya,m,p)
%---------------------------------------------------------------------------

  % RK4   Runge-Kutta solution for y' = f(t,y) with y(a) = ya .
    % Sample call
    %   [T,Y] = rk4('f',a,b,ya,m)
    % Inputs
    %   f    name of the function
    %   a    left  endpoint of [a,b]
    %   b    right endpoint of [a,b]
    
    %   ya   initial value
    %   m    number of steps
    %   p    Kuramoto function arguments
    % Return
    %   T    solution: vector of abscissas
    %   Y    solution: vector of ordinates
    %
    % NUMERICAL METHODS: MATLAB Programs, (c) John H . Mathews 1995
    % To accompany the text:
    % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
    % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U . S . A .
    % Prentice Hall, Inc .; USA, Canada, Mexico ISBN 0-13-624990-6
    % Prentice Hall, International Editions:   ISBN 0-13-625047-5
    % This free software is compliments of the author .
    % E-mail address:      in % "[e-mail]"
    %
    % Algorithm 9.4 (Runge-Kutta Method of Order 4) .
    
    % Section	9.5, Runge-Kutta Methods, Page 460
    %---------------------------------------------------------------------------
        h = (b - a)/m;
        T = zeros(1,m+1);
        size(ya,1)
        Y = zeros(size(ya,1),m+1);
        T(1) = a;
        Y(:,1) = ya;
        for j=1:m
          tj = T(j);
          yj = Y(:,j);
          k1 = h*feval(f,tj,yj,p);
          k2 = h*feval(f,tj+h/2,yj+k1/2,p);
          k3 = h*feval(f,tj+h/2,yj+k2/2,p);
          k4 = h*feval(f,tj+h,yj+k3,p);
          Y(:,j+1) = yj + (k1 + 2*k2 + 2*k3 + k4)/6;
          T(j+1) = a + h*j;
        end
end

[Walter Roberson]

Are you certain that you want to use the matrix-right-divide operator, / there ? A/B is similar to A*pinv(B) but with some additional restrictions, and effectively does a least-squared fit. Are you sure that is what you want to do?

Or do you perhaps want to use the element-by-element division operator, ./ ?

[H_W]

@Walter Roberson I've tried A*pinv(B), but obtained non-real sol, and after several iterations it returns NaN + NaNi.

Answer 1

[Bruno Luong]

Replace ^2 by .^2 (with the dot) in fe if you want to square element wise.

There migjt be some other operation like * or / that you need to check if that is correct (I suspect not).