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

5 views (last 30 days)
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 % "mathews@fullerton.edu"
%
% 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
  2 Comments
H_W
H_W on 4 Nov 2022
Edited: H_W on 4 Nov 2022
@Walter Roberson I've tried A*pinv(B), but obtained non-real sol, and after several iterations it returns NaN + NaNi.

Sign in to comment.

Answers (1)

Bruno Luong
Bruno Luong on 4 Nov 2022
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).

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!