ODE45 generates undesired matrix of NaN entries

I am trying to use ODE45 to solve a system of 6 differential equations with 6 unknowns. I have defined them in a function but when I try to integrate using ODE45 I get a matrix of all NaN entries except for the very first row. I know that NaN comes from something like dividing by 0 but I dont see where that could possibly occur in my code. Here is my code so far.

% Function for 6 diff eq's
function int = brk(~,x3)
a3 = 0.140563781397801;
b3 = -2.931811125246462;
uk = 0.089552917882572;
Mp = .1; Md = 10; C = 4.63e10; kp = 39500; kd = 3.95e8; vi = .01;
F = 1.2;
% 6 equations
int(1) = x3(4);
int(2) = -kp*x3(1)/Mp + (C/Mp)*((x3(2)-x3(3))^(3/2))*(uk+a3*exp(b3*(vi-x3(4))));
int(3) = x3(5);
int(4) = F/Mp - (C/Mp)*((x3(2)-x3(3))^(3/2));
int(5) = x3(6);
int(6) = (C/Md)*((x3(2)-x3(3))^(3/2)) - kd*x3(3)/Md;
int = int(:);
end

And here is the script where I try to use ODE45 to solve the function

% Integration for position differential equations
clear all
a3 = 0.140563781397801;
b3 = -2.931811125246462;
uk = 0.089552917882572;
Mp = .1; Md = 10; C = 4.63e10; kp = 39500; kd = 3.95e8; vi = .01;
F = 1.2;
h3 = 0.01;
% Initial conditions
icond = [0 0 ((F/C)^(3/2))+F/kd 0 F/kd 0 ];
tstep = 0:h3:.5;
[t3,ode] = ode45('brk',tstep,icond);
% Plot results for x1
plot(t3,ode(:,1))

If anyone could help me figure this out that would be amazing. I have been trying to fix it forever now and cant get it to work.

Answer 1

[Jan]

Use the debugger to stop, when a NaN occurs:

dbstop if naninf

Type this in the command window. Now run the code again until it stops at the first occurring NaN. Now you can check the current values.

By the way, in older Matlab versions clear all cleared the debugger status also. Everything which impedes debugging is an enemy of the programmer. Therefore I strongly recommend to omit this brute clearing, which wastes time only. Prefer to use functions, if you want to keep the workspace clean.

Defining the function to be integrated by a char vector is outdated for 15 years now, but supported for backward compatibility. Use a function handle instead:

ode45(@brk, ...

[james poston]

I got rid of the clear all at the beginning and used the debugger stop and it opened up a new script which stopped on this line:

 case 'getvariableclass'
        throwError = true;
        try
            if nargin > 1 && isa(varargin{1},'com.mathworks.mlservices.WorkspaceVariable') && length(varargin{1}) == 1
                varargout = {class(currentValue)};
            elseif nargin > 1
                % try to eval the variable here, instead of when calling
                % the function, so that we can catch any errors and prevent
                % error beeping
                v = evalin('caller', varargin{1}); % <- This line is where it stopped
                varargout = {class(v)};
            else
                varargout = {''};
            end

However I am not sure what this means or how to change my code. What does this debug info mean and how do I use it?

[Jan]

@james: Ups? Which function is this? Actually I'd expect the debugger to stop inside brk(). From whereis this function called?

If I run your code with dbstop if naninf, it stops at:

int(2) = -kp*x3(1)/Mp + (C/Mp)*((x3(2)-x3(3))^(3/2))*(uk+a3*exp(b3*(vi-x3(4))));

when x3(4) equals 6.3219e+12 + 4.2089e+12i. Then exp term replies -Inf + inf i.

I've added the variable t in the function call:

function int = brk(t,x3)

and find, that the Inf value occurs at t=1.6746e-05, and this seems to be the first step of the integration already. Therefore I assume the formula is wrong or has a pole at the initial value.

Answer 2

[James Tursa]

This is often an initial condition issue with the specific DE's involved. In your code, you have

icond(2) = 0
icond(3) = ((F/C)^(3/2))+F/kd = 3.0380e-09

Then in your derivative function, you have this expression:

(x3(2)-x3(3))^(3/2)

So for the very first call to your derivative function, that's going to result in a negative value being raised to the 1.5 power ... a complex result. I am pretty sure this is not what you want.

You need to go back and look at your original DE's (which you didn't provide to us) and your initial conditions again.