RuntimeWarning: invalid value encountered in double_scalars app.launch_new_instance()

Question

I am applying Euler's method to solve a differential equation. This is my code:

def f(x, y):
    return ((x**(2))*y)/((x**(4)) + (y**(4)))

di = 0.01
I = 100

x = np.linspace(-I, I, int(I/di) + 1)
w = np.zeros(len(x))

x[0], w[0]
for i in range(1, len(w)):
    w[i] = w[i - 1] + f(x[i - 1], w[i - 1])*di

plt.plot(x, w, label='approximation' )
plt.xlabel("x") 
plt.ylabel("y")
plt.show()

When I run the code I have this Warning:

"C:\Users\USER\Anaconda3\lib\site-packages\ipykernel__main__.py:3: RuntimeWarning: invalid value encountered in double_scalars app.launch_new_instance()"

I would like to now how to solve it and make it work.

Answer 1

Your code is running into Divide by Zero Error. Try this to convince yourself:

>>> def f(x, y):
...     return ((x**(2))*y)/((x**(4))+(y**(4)))
...
>>> I, di = 100, 0.01
>>> x = np.linspace(-I, I, int(I/di) + 1)
>>> w = np.zeros(len(x))
>>> i = len(x)//2 + 1
>>> i
5001
>>> f(x[i-1], w[i-1])
nan

It clearly emerges from the interactive session above that when i takes the value 5001 in the for loop, f(x[i-1], w[i-1]) yields nan. There are several solutions to this issue. For instance, in order to avoid NaN values you could check whether the denominator of the fraction returned by f() is zero prior to perform the division. If it is, you should return a conventional value of your choice (for example 0) instead of the result of the division. The following snippet implements such approach though a conditional expression:

def f(x, y):
    return (0 if x==0 and y==0 else float(x**2*y)/(x**4 + y**4))

Alternatively, you could disable run time warnings (but if you do so you need to be aware of the potential risks) by including this code in your script:

import warnings

def f(x, y):
    with warnings.catch_warnings():
        warnings.simplefilter('ignore')        
        return ((x**(2))*y)/((x**(4)) + (y**(4)))

---

The proposed workarounds avoid the RuntimeWarning but don't get your code working as you expect. Indeed, the calculated solution w is a vector in which all the elements are zero. I guess the reason for your code not being working properly is that you missed to assign w[0] an initial value different to 0.

For example, if you simply add this line just before the for loop:

w[0] = 0.5

you get this (apparently meaningful) curve rather than a flat plot.

Hope this helps!