You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
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?
25 views (last 30 days)
Show older comments
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:
Thanks a lot!
13 Comments
dpb
on 3 Jul 2020
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...
dpb
on 3 Jul 2020
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
on 3 Jul 2020
Edited: Sim
on 6 Jul 2020
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
on 3 Jul 2020
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
on 3 Jul 2020
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
on 3 Jul 2020
I think that's what happens to socks that get lost in the laundry ... they morph into coat hangers.
Sim
on 6 Jul 2020
Edited: Sim
on 6 Jul 2020
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
on 6 Jul 2020
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.
John D'Errico
on 7 Jul 2020
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
on 10 Jul 2020
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)
Answers (0)
See Also
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)