2

I'm trying to do a bifurcation analysis of a second order dynamical system using MatCont package of MATLAB. I'm fairly new to MatCont and ODE solvers in general.

The system is the following (as written in MatCont):

x1' = R*x1/L + x2/L
x2' = a+x2/D - (b*x2^3)/D - x1/D

Where R,L,D,a,b are parameters. When I start the integration focusing on just one single variable (my aim is to plot one variable against the bifurcation parameter R), I get the following error:

"current step too small".

I've tried different values for InitStepSize, MinStepSize and MaxStepSize but the error doesn't go away, so I'm not able to perform the analysis.

The solver I'm using is ode45. Any idea how to fix this issue??

More details:

  1. First 3 derivatives are computed symbolically
  2. Initial point: (t=0, x1=0, x2=0, R=27, L=1.2e-08, D= 1.2e-11, a=0.04, b=0.18)

This is the output I get from matlab:

enter image description here

first point found
tangent vector to first point found
Current step size too small (point 1)
elapsed time  = 2.2 secs
npoints curve = 1
first point found
tangent vector to first point found
Current step size too small (point 1)
elapsed time  = 0.4 secs
npoints curve = 1

If you need more details, please ask.

ouflak
  • 2,458
  • 10
  • 44
  • 49
AlePalu
  • 73
  • 1
  • 8
  • Does the error originate from inside MatCont or ode45? There are two stepping where the step size could get too small to further advance the variable, the R axis and the time axis. // What part are the options for, MatCont or ode45? // It would be helpful if you could post more code and the full error message. – Lutz Lehmann Dec 07 '20 at 12:09
  • @LutzLehmann thank you for your answer. I've edited the question with the output I get from MatCont. Unfortunately I've no code to show, simply because I just use the graphical interface offered by MatCont, I've just specified the model and set the parameters. As for the choice of the solver, ode45 is an option available from the MatCont interface, so I've simply selected it. The options I've reported in the question are hence for MatCont, since in no way I've interacted directly with the solver. Can be the choice of the solver a problem? I mean, changing the solver could make some difference? – AlePalu Dec 07 '20 at 15:30
  • Yes, you could use one of the implicit solvers, in some cases of stiff systems they allow larger step sizes. But your output points more to a root finder with line-search. // You could try to scale the time by a factor `1e9` so that the constants `L,D` get values on a more appreciable scale, `L=12` and `D=0.012`. Sometimes algorithms are designed with heuristics that are calibrated for simple scales and do not contain measures to deal with problems with widely different scales. – Lutz Lehmann Dec 07 '20 at 16:23
  • @LutzLehmann thanks to your advice I was able to solve the issue! :) Just scale the parameters as suggested by you and the solver (still ode45) is now able to perform the integration correctly. If you want you can answer to my question so that your suggestion will be helpfull to other people (there is little documentation on this problem on internet, expetially how to solve this issue). – AlePalu Dec 08 '20 at 10:39

1 Answers1

2

When implementing an algorithm, one has to decide how general one wants to make it. One of the most intuitive assumption that can go wrong when confronted with "real-world" problems is that all variables have about the same weight, that all regions of interest are non-flat, for instance nearly circular or spherical in the unmodified Euclidean norm or scalar product.

There are two main ways to achieve this assumption, pass a scale vector that is used to construct an adapted norm and scalar product, or just pass such a scalar product, or just demand that the user already has scaled their problem in such a way that the range of the values of each variable is about 0.01 to 1000 or so.

A sign that this might not be the case are extreme values in "connection" coefficients in the equations like here L and D. While ode45 might still gracefully deal with this, the non-linear root finder behind the bifurcation analysis may not. With an easy change in the time scale by a factor of 1e9 (that is, instead of seconds as time unit you go to nano-seconds), the time scale becomes comparable to the state space scales. In code this change is achieved by simply multiplying L,D with this factor, giving L=12 and D=0.012. As reported, this actually already was sufficient to make the analysis program work.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51