0

I am trying to solve a system of non linear differential equations by the RK4 method. One of my parameter in that equations is a very small number on the order of 10^{-121} thus to work with more precision I am using Decimal package. But I am stuck with the "decimal.InvalidOperation" warning.

def V(x,y,N):
v = np.divide(y**2, (1-x**2-y**2)) * 3* H02 * (Om1 * np.exp(-3*N) + Om2 * np.exp(-4*N))
return v

def f1 (x,y,l,q):
  return -3*x + l * np.sqrt(3/2) * y**2+ 3/2 * x * ( 2*x**2 + q*(1 - x**2 y**2)) 

def f2 (x,y,l,q): 
  return - l * np.sqrt(3/2) * y * x + 3/2 * y * (2 * x**2 + q*(1 - x**2 - y**2))

def f3 (N, x, y, l, g, M, al): 
  exp1,   G,   m  = decimal.Decimal(-2*g).exp(),    decimal.Decimal(g),         decimal.Decimal(M)
  xx,  a,   ll    = decimal.Decimal(x), decimal.Decimal(al),    decimal.Decimal(l)
  VV              = decimal.Decimal(V(x,y,N)) 
  log             = decimal.Decimal(VV/(m*exp1) + 1).ln()
  Tanh            = 1/G * log  - 1
  # gamm            = -(G*(1 - Tanh**2) - 2 * Tanh)/(np.sqrt(6*a)*ll)
  return  float(ll *xx/np.sqrt(a)* (G*(1 - Tanh**2) - 2 * Tanh + np.sqrt(6*a)*ll))

def R_K_4(xi,yi,li, g, M, al):
  N = int(round((lnaf-lnai)/dlna))
  lna = np.linspace(lnai, lnaf, N)
  x = np.zeros(N) 
  y = np.zeros(N)
  l = np.zeros(N)
  x[0],y[0],l[0] = xi,yi,li
  for i in range(0,N-1):

    q    = 1 + 1/(3*(1 + 5780*np.exp(lna[i])))

    kx1 =  f1( x[i], y[i], l[i], q )
    ky1 =  f2( x[i], y[i], l[i], q )
    kl1 =  f3( lna[i], x[i], y[i], l[i], g, M, al)
    
 
    kx2 =  f1(                     x[i]+kx1*dlna/2,  y[i]+ky1*dlna/2,  l[i]+kl1*dlna/2, q )
    ky2 =  f2(                     x[i]+kx1*dlna/2,  y[i]+ky1*dlna/2,  l[i]+kl1*dlna/2, q )
    kl2 =  f3( lna[i] +  dlna/2,   x[i]+kx1*dlna/2,  y[i]+ky1*dlna/2,  l[i]+kl1*dlna/2,       g, M, al )
    

    kx3 =  f1(                     x[i]+kx2*dlna/2,  y[i]+ky2*dlna/2,  l[i]+kl2*dlna/2, q )
    ky3 =  f2(                     x[i]+kx2*dlna/2,  y[i]+ky2*dlna/2,  l[i]+kl2*dlna/2, q )
    kl3 =  f3( lna[i] +  dlna/2,   x[i]+kx2*dlna/2,  y[i]+ky2*dlna/2,  l[i]+kl2*dlna/2,        g, M, al )

    kx4 =  f1(                     x[i]+kx3*dlna,    y[i]+ky3*dlna,    l[i]+kl3*dlna,   q )
    ky4 =  f2(                     x[i]+kx3*dlna,    y[i]+ky3*dlna,    l[i]+kl3*dlna,   q )
    kl4 =  f3( lna[i] +  dlna,     x[i]+kx3*dlna,    y[i]+ky3*dlna,    l[i]+kl3*dlna,          g, M, al )

    x[i+1]  = x[i] +dlna/6*(kx1 + 2*kx2 + 2*kx3 + kx4)
    y[i+1]  = y[i] +dlna/6*(ky1 + 2*ky2 + 2*ky3 + ky4)
    l[i+1]  = l[i] +dlna/6*(kl1 + 2*kl2 + 2*kl3 + kl4)

  return x,y,l

and when I call the solver as

H02 = 1.44*1e-122
Om1, Om2 = 0.31, 5.45*1e-5
lnai,  lnaf, dlna  = -15, 10, 1e-2
x11, y11, l11 = R_K_4(0, 2.4e-11, -0.919, 128, 1e-10, 7/3)

I am getting the above mentioned error. But I think the error takes place before that step but I am unable to identify where exactly something went wrong. Please help me with this frustration!

user280016
  • 19
  • 2
  • Can you post the full error traceback? – MegaIng Apr 30 '23 at 12:15
  • Also, please provide an actual [MRE]. The code you currently have produces `NameError`s because of the undefined name `lnaf`. – MegaIng Apr 30 '23 at 12:18
  • Thanks for pointing it out! Now, I hope it is reproducible. – user280016 Apr 30 '23 at 13:37
  • The problem is that your potential function returns a negative value. `ln(-x)` is not defined, neither for decimal nor for floats. – MegaIng Apr 30 '23 at 14:04
  • If you want to use decimal, everything should be decimal based: no `numpy`, no raw floats anywhere in the source code – MegaIng Apr 30 '23 at 14:08
  • I noticed that myself too but as I was saying the potential becoming negative happens during the iteration in my R_K_4(), so something goes wrong during the evaluation before it comes to evaluate my potential. Otherwise, the code is expected for 100% to produce a sensible result. – user280016 Apr 30 '23 at 14:11
  • to your last comment, I solved this problem with 2D system instead of 3D (as it is now) and there I used decimal only to the functions that are sensitive to small changes (not to everything since it slowed down the computation time too much); it worked just fine but when I want to repeat the same for my 3D equations I have this error message. – user280016 Apr 30 '23 at 14:16

0 Answers0