I have to make a program using Euler’s method for the “ball in a spring” model
from pylab import* from math import* m=0.1 Lo=1 tt=30 k=200 t=20 g=9.81 dt=0.01 n=int((ceil(t/dt))) km=k/m r0=[-5,5*sqrt(3)] v0=[-5,5*sqrt(3)] a=zeros((n,2)) r=zeros((n,2)) v=zeros((n,2)) t=zeros((n,2)) r[1,:]=r0 v[1,:]=v0 for i in range(n-1): rr=dot(r[i,:],r[i,:])**0.5 a=-g+km*cos(tt)*(rr-L0)*r[i,:]/rr v[i+1,:]=v[i,:]+a*dt r[i+1,:]=r[i,:]+v[i+1,:]*dt t[i+1]=t[i]+dt #print norm(r[i,:]) plot(r[:,0],r[:,1]) xlim(-100,100) ylim(-100,100) xlabel('x [m]') ylabel('y [m]') show()
I keep getting this error:
a=-g+km*cos(tt)*(rr-L0)*r[i,:]/rr RuntimeWarning: invalid value encountered in divide
I can’t figure it out, what is wrong with the code?
I think your code is trying to “divide by zero” or “divide by NaN”. If you are aware of that and don’t want it to bother you, then you can try:
import numpy as np np.seterr(divide='ignore', invalid='ignore')
For more details see:
Python indexing starts at 0 (rather than 1), so your assignment “r[1,:] = r0” defines the second (i.e. index 1) element of r and leaves the first (index 0) element as a pair of zeros. The first value of i in your for loop is 0, so rr gets the square root of the dot product of the first entry in r with itself (which is 0), and the division by rr in the subsequent line throws the error.
To prevent division by zero you could pre-initialize the output ‘out’ where the div0 error happens, eg
np.where does not cut it since the complete line is evaluated regardless of condition.
example with pre-initialization:
a = np.arange(10).reshape(2,5) a[1,3] = 0 print(a) #[[0 1 2 3 4], [5 6 7 0 9]] a/a # errors at 3/0 out = np.ones( (5) ) #preinit np.divide(a,a, out=out, where=a!=0) #only divide nonzeros else 1
You are dividing by
rr which may be 0.0. Check if
rr is zero and do something reasonable other than using it in the denominator.