[SciPy-Dev] Choosing the result of dtype for integration

[yotam vaknin]

Thanks David for your response!

I understand your desire to keep the surprising parts of the complex domain (sqrt, log) on an opt-in basis. But a lot of the non problematic parts are not opt-in at the moment. Functions with real images (with respect to the reals, like exp, sin, erf) don’t have a complex opt-in functionality. I would expect integration to fall into this category, since keeping the status quo results in behavior which is very unintuitive (as oppose to log(-1)). 

If I understand everything correctly, the current implementation of the sqrt function shouldn’t result in any problem, since sqrt doesn’t return complex values when one doesn’t ask for them, keeping the complex functionality on an opt-in basis. 

Regards,
Yotam Vaknin

> On 28 Feb 2020, at 3:46, David Hagen <david at drhagen.com> wrote:
> 
> 
> I am all for improving handling of complex integration, but I think the current behavior matches how SciPy handles real and complex domains, namely by making the user opt into the complex domain. In Matlab log(-1.0) is 3.14i. In NumPy, log(-1.0) is NaN. NumPy and SciPy require the user to be explicit about working in the complex domain. Matlab happily takes you from the real domain to the complex domain without warning. I have done a lot of ODE integration in Matlab, and I always disliked how easily it accepted complex numbers. If I had a log somewhere in my RHS and it accidentally got a negative number, the whole thing would keep running and just return a giant block of complex nonsense. Now, this is a bit different because it is a lot harder to accidentally get a complex RHS, which is why I am only weakly in favor of maintaining the status quo.
> 
>> On Thu, Feb 27, 2020 at 5:45 AM yotam vaknin <tomirendo at gmail.com> wrote:
>> Hello list,
>> 
>> I hope I'm now visiting a very old topic, but I wanted to suggest a simple change to the integration mechanism.
>> 
>> At this moment, the dtype of the result of an integration is the dtype of the initial condition. For example, when using ״solve_ivp" to solve the initial value problem 
>> y0=np.array([1,0] )
>> dy/dt= -i H at y0
>> which is a very a standard calculation in Quantum Mechanics, the integrator will return nonsense, since the process involves complex numbers, and the initial conditions were (wrongly assigned to be) real. The solution is of course y0=np.array([1,0], dtype=complex), but this is very hard to figure out to new users coming from, for e.g., Matlab.
>> 
>> I wanted to suggest that instead of choosing the dtype using the initial condition, that the resulting dtype will be the dtype of the derivative (or better yet, derivative*dt + initial condition). This will solve the problem, and introduce minimal confusing.
>> 
>> Thanks for all of your great work!
>> 
>> Regards,
>> Yotam Vaknin
>> 
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev at python.org
>> https://mail.python.org/mailman/listinfo/scipy-dev
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at python.org
> https://mail.python.org/mailman/listinfo/scipy-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20200228/b1496925/attachment.html>