0

I compute something in python with sympy but this require too much time. The code looks like this:

def compute_aux_lyupanov_val_k(y, A, B, C, D, k):
  res = []
  res.append(y)
  for i in range(k-1):  
    y = (C + D * y) / (A + B * y)
    res.append(y)
  return res


def compute_lyupanov_val(y, A, B, C, D, n):
  res = 0
  y_values = compute_aux_lyupanov_val_k(y, A, B, C, D, n) # values for y'_0. y'_1...

  for i in range(n):
    T = ( (A + B * y_values[i])**2 + (C + D * y_values[i])**2  ) / (1 + y_values[i]**2)
    res = res + math.log(T, 10) # default is base e = 2.71

  res = res * (1/(2*n))
  return res

n = 1000
y_inital = 0.01

beta_set_values = np.linspace(0, 6, 1000)
alpha_set_values = np.linspace(0, 5, 1000)
points = it.product(beta_set_values, alpha_set_values)

for point in points:
  alpha_val = point[1]
  beta_val = point[0]

  A_val = A.subs([  (beta, beta_val), (alpha, alpha_val) ])
  B_val = B.subs([  (beta, beta_val), (alpha, alpha_val) ])
  C_val = C.subs([  (beta, beta_val), (alpha, alpha_val) ])
  D_val = D.subs([  (beta, beta_val), (alpha, alpha_val) ])

  lyupanov_val = compute_lyupanov_val(y_inital, A_val, B_val, C_val, D_val, n)

  if lyupanov_val > 0:
    plt.scatter([beta_val], [alpha_val], color = 'r', s=1)


plt.show()

Maybe it seems long, but let me explain clearly: A, B, C, D are expression which can be evaluated depending on alpha, beta values. For a set of points (alphas and betas) I compute a value (lyupanov_val) and scatter the (beta, alpha) point only if computed value lyupanov_val is greater than 0. The computations which happened in auxiliary functions (compute_aux_lyupanov_val_k and compute_lyupanov_val) are only numerical (A, B, C, D and other parameters are scalars now - inside of these functions, not symbols)

I think the expensive computation part happen because of the substitution at every iteration. I know that sympy provide a special plot function which can speed up the computations sp.plot_implicit(expresion > 0, (beta, 0, 6), (alpha, 0, 5) but the problem with this approach is the fact that I don't use a direct expression (e.g alpha + beta / 2), I use auxiliary functions for numerical computations (e.g compute_aux_lyupanov_val_k)

If I used

sp.plot_implicit(compute_lyupanov_val(y_inital, A, B, C, D, n) > 0, (beta, 0, 6), (alpha, 0, 5))

then the python raise TypeError: can't convert expression to float

Any ideas?

dacian
  • 95
  • 1
  • 8
  • 2
    Calling `plt.scatter` in a loop is likely to be the most expensive part. Collect the points in a list and call `plt.scatter` once at the end. – Oscar Benjamin Jan 08 '23 at 12:39
  • If I read this correctly, `A,B,C,D` are sympy expressions, but `A_val, etc` are scalars. And `compute_lyupan` is plain python, just numeric. I wonder if you can `lambdify` `A,B,C,D`, and given them the arrays produced by `np.meshgrid` (or `mgrid`) from your `linspace` arrays. And rewrite `compute_lyupan` to work with arrays. The only `>` that would prevent numpy "vectorization" is in the choice of values to plot. – hpaulj Jan 08 '23 at 17:19
  • The iteration in `compute_aux_lyupanov_val_k` might be hard to "vectorize", since at each step `y` is a function of the previous `y`. I wonder if it can be changed to use something like `np.cumprod`. – hpaulj Jan 08 '23 at 17:24

1 Answers1

0

subs does a symbolic substitution but the result is still a symbolic expression. May be what you want is instead

A_val = A.subs([  (beta, beta_val), (alpha, alpha_val) ]).evalf()

and so on for the 4 variables.

Even better would be to use

A_val = A.evalf(subs={beta: beta_val, alpha: alpha_val})

(subs does a replacement and a symbolic simplification; if what you want is instead just the numeric value a direct computation is better).

6502
  • 112,025
  • 15
  • 165
  • 265
  • 1
    yes, I want only numerical value; seems like A.evalf(...) convert to *numerical sympy* (the type of the resulted value is also related to sympy), so I need to use also float (A.evalf(...)) but thank you for suggestion – dacian Jan 08 '23 at 12:51