Getting complex solutions with MATLAB's ode45

Question

I'm trying to solve the ODE below using Matlab's ode45:

$$\frac{dy}{dx}=\sqrt{A\cdot \sin(y) + B\cdot \cos(y)}$$

where $0\le x \le 1$.

The Matlab code is:

A = -2.45;
B = 2.50;
tspan  = 0:0.01:1;
odefun = @(t,y) sqrt(A*sin(y) + B*cos(y));
[t,y] = ode45(odefun, [0 1], 0);

The constants $A$ and $B$ are chosen such that inside the square root it is always non-negative. The solution seems fine for most values of $x$, but for the last few points close to $x=1$ I get complex $y$ with very small imaginary parts. I then took the real or the absolute value of $y$, but then I would get negative value inside the square root.

Should I try other ODE solvers? Or is there anything else I could do to eliminate the complex solutions?

Answer 1

This first-order ODE is such that $dy/dx \geqslant 0$. Therefore, a solution $y$ must be an increasing function of $x$. Here, we start with the initial condition $y(0)=0$, where the derivative satisfies $dy/dx \approx 1.6 > 0$. The derivative $dy/dx$ becomes complex when its square becomes negative. For this to happen, an equilibrium value (i.e. where the derivative vanishes) must have been reached. Here, the smallest positive equilibrium value is \begin{equation} y^* = {-\arctan(B/A)} \approx 0.796 \, . \end{equation} This can be observed when using Matlab: the value $y^*$ is reached numerically just before $x=1$. The choice $A=2.45$ and $B=2.5$ does not avoid $dy/dx$ to become complex.

Note: the Picard-Lindelöf theorem applies here. If $x^*$ is such that $y(x^*)=y^*$, we know that the solution $y$ over $[0, x^*[$ is unique. Thus, it can be approximated numerically with classical methods such as ode45.

Answer 2

In reviewing my previous answer it became apparent that the solution there is not the most simple one.

The most straight-forward solution is to extend the domain of the differential equation without changing it on the given domain, so that reaching and slightly crossing the boundary of the given domain does not lead to non-real values for the derivatives and then the values of the numerical solution. This in the main part is achievable here by replacing the square root with an extension to the whole real line. For the solution to stabilize at the boundary one can take the odd extension $$ r(u)={\rm sign}(u)\sqrt{|u|}. $$ Numerical solvers find large or infinite values for derivatives difficult to deal with, leading to very small step sizes. To mollify this one can replace some small segment around $u=0$ with a line like in $$ r(u)=\frac{u}{\max(\epsilon,\sqrt{|u|})} $$ or similar constructs. To get large step sizes in the asymptotic segment of the solution one will need implicit solvers. Using the implicit "Radau" solver gives the following picture where the internal steps are noted with crosses: