Is it possible that ode23 or ode45 could result in complex numbers instead of real ones? If Yes, what are the possible reasons for getting complex numbers once you would expect real ones?

Hi, I perfromed ode23 for solving a differential equation, but after some timesteps I started to get complex numbers (i.e. the imaginary part) as result, instead of real numbers. Based on my equations/calculations, I should always get real numbers, and I do not understand if/why those solvers could give imaginary numbers.

Therefore, Is it possible that ode23 or ode45 could result in complex numbers instead of real ones?

If Yes, what are the possible reasons for getting complex numbers once you would expect real ones?

More info on the ODE here:

https://it.mathworks.com/matlabcentral/answers/558265-error-in-using-ode23t-to-solve-a-differential-equation-on-a-graph-network?s_tid=mlc_ans_email_ques

Thanks a lot!

[dpb]

Having square roots or exponentials or any number of other functional forms in the system of equations could lead to the symptom if.

We have no way to know what your ODE looks like...

[Sim]

Hi, thanks for your answer :) ...I forgot to post the ODE.... that 's in my previous question here:

https://it.mathworks.com/matlabcentral/answers/558265-error-in-using-ode23t-to-solve-a-differential-equation-on-a-graph-network?s_tid=mlc_ans_email_ques

[dpb]

If it's the same Q?, don't start new thread; if it's different Q? but happens to use same ODE, then put the info here...

[Sim]

I guess "Q = question".... It is a different question but with the same ODE:

% network
s = [1 1 2 2 3 3 4 5 5 6  6 7  7  8  9  9 10 10 11 11 12 13 14 15];
t = [2 5 3 6 4 7 8 6 9 7 10 8 11 12 10 13 11 14 12 15 16 14 15 16];
G = graph(s,t);
N = numnodes(G); % 16 nodes
e = table2array(G.Edges);
n = table2array(G.Nodes);
% networks coordinates (x,y,z)
x = ...
y = ...
z = ...
G.Nodes.X = x'; G.Nodes.Y = y'; G.Nodes.Z = z';
% inputs/initial conditions/parameters
initialnode = 1;
h = repelem(0.01,N); 
h(initialnode) = 0.5; 
w = repelem(1,N); 
Q(initialnode) = 0; 
lambda = 0.3;
% ode23t solver
tspan = [0 3600]; % t = time
opts = odeset('MaxStep',3600);
[T,H] =  ode23t(@(t,h) MYODE(t,h,G,N,e,n,z,initialnode,w,Q,lambda),tspan,h,opts);
% MYODE 
function dhdt = MYODE(t,h,G,N,e,n,z,initialnode,w,Q,lambda),tspan,h,opts);
    for i = 1 : N
        j = randi([1 N]);
        % some operation/equation (I use cells for Q, since it could "store" several numbers)
        Qcell{i} = ... % some operation/equation
        dhdt_tmp{i} = cell2mat(Qcell{i}) / function_ij; % differential equation
    end
    dhdt = cell2mat(dhdt_tmp);
end

[dpb]

Qcell{i} = ... % some operation/equation

Well, we still have no idea what might be in Qcell or how it was calculated. Wherever/whatever it is, there's where the complex variables will be coming from.

[John D'Errico]

I could swear I recall seeing someone ask a question recently about an ODE that was producing complex results. The ODE had terms in it like x^n, where n was non-integer. Note that if x ever goes negative, you will generate complex results. But your variable might have that happen in an intermediate result. Once complex enters te computations, complex numbers tend to propagate. Complex numbers (and NaNs) are like wire coat hangers. They just breed like crazy.

[James Tursa]

I think that's what happens to socks that get lost in the laundry ... they morph into coat hangers.

[Sim]

Thanks to everyone for your insightful answers :-)

@dpb, this is Qcell{i}

% c1, c2, c3 are constant numbers and h is the unknown quantity of the differential equation
A = (c1 * h) / (c2 + h);
B = sqrt(c3 * A);
Qcell{i} = B * h;

@John D'Errico, Yes, in my case, the unknown quantity "h" goes suddenly negative and after that a moltitude of complex numbers occurs as in a chain reaction!

So, probably, in ODE23, my unknown quantity goes negative here (Ordinary Differential Equation Solvers ODE23 and ODE45)

% https://blogs.mathworks.com/cleve/2014/05/26/ordinary-differential-equation-solvers-ode23-and-ode45/
yn+1 = yn + h/9 * (2*s1 + 3*s2 + 4*s3)

??!

@James Tursa, Hahaha, its true! :)

[John D'Errico]

Ah. I was right. Those blasted complex numbers. You see? They do breed, making little baby complex numbers at first. Should I check my sock drawer though? Perhaps the dust bunnies I found in there were really little baby socks?

Seriously, how do you fix it? You could most simply run the ODE solver until the point where the variable of interest goes negative. Use the odeevents option to identify that point. Then stop the solver. Restart it again, now going in the correct direction.

[Sim]

Yes, I'll try odeevents, many thanks :) ....Also, I would like to accept all your answers, but I cannot see the "Accept" button anymore...

[Steven Lord]

The NonNegative option may also be of use. See this documentation page for more information.

[John D'Errico]

Sometimes when we don't know what the answer is, a group of people end up just making comments. Along the way, sometimes we might end up helping you, finding a solution as if by committee. If that worked, then great, you have a solution, and all is good in the world. What matters in the end is you have a solution to the problem.

Just watch out for those multiplying socks and hangars. :)

[Sim]

Thanks a lot to you all :)

@Steven Lord, thanks for your suggestion :) However, I do not see a so much difference by using the NonNegative option...

@John D'Errico, Yes you all gave me the right advices, thanks a lot :)

I think I partially fixed the issue of complex numbers by putting an abs() function here (please see my previous comment to @dpb, about the Qcell{i})

% c1, c2, c3 are constant numbers and h is the unknown quantity of the differential equation
A = abs((c1 * h)) / (c2 + h);
B = sqrt(c3 * A);
Qcell{i} = B * h;

However, once I run the code - and by using ode15s instead of the previous ode23t, since a bit faster than ode23t - I get this plot, which represents my unknown quantity "h" versus time "t", for different parts of my network:

At the beginnning, up to about t=13, everything looks correct to me, i.e. positive smooth curves (curve = equation's solutions).. However, after that time, some curves go negative and some get a spike..... Why?! :)

(Please let me know if I should open a new thread about this last question)