0

I have an implicit equation, like this:

(a1 X + b1 Y + m)*(a2 X + b2 Y + m)*(a3 X + b3 Y + m) - c = 0,

a1,a2,b1,b2,a3,b3 are certain value, c is a variant. According different c, I need to solve it to get a set of (x,y), which I will use to integrate.

The listed function I have in practice is much more complex than this, so I am confused as to why when I put this equation into the website desmos to draw this implicit function, I can get the solutions that satisfies this function and I would like to know why this is so fast for desmos and then if there is a better way to find these solutions

I using polar coordinate to solve this problem

c_range = np.linspace(0,c_max, 1000)
theta_range = np.linspace(0,pi, 1000)
for i in range(1000):
    if i ==0:
        c = c_range[i]
        for j in range(1000):
            theta = theta_range[j]
            r = fsolve(@func, r0, args=(c, theta))
            radius[i][j] = r
            r0 = r
    else:
        c = c_range[i]
        r0 = radius[i-1]
        r = fsolve(@func, r0, args=(c, theta_range))
        radius[i]= r
dali
  • 1
  • 1
  • The solution of an implicit equation is a contour in the x-y plane. There exist simple heuristics to follow such a contour on a grid. So the question is, what do you do that is so slow? – Lutz Lehmann May 17 '22 at 09:34
  • Yes it is a series of contour lines according different c, I need these points (x,y) on the contour lines, to integrate over x. I use for loop to get these coordinates, the similar code I have put up. what is heuristics? sorry I am new to python and stack overflow – dali May 17 '22 at 16:23
  • This looks already quite optimal, in the sense of general ideas. You could improve the next guess via extrapolation. Using fsolve has a large overhead, using it vectorized is perhaps not a good idea as the full Jacobian gets computed where you know that it is always diagonal. Also, the overhead of fsolve might be rather large, this could perhaps be improved with the scalar root-finder `find_root`. // The keyword I was missing above is [marching squares](https://en.wikipedia.org/wiki/Marching_squares) – Lutz Lehmann May 17 '22 at 17:09
  • Yes it costs a lot of time. But do you mean using scipy.optimize.root? I have tried this but it costs more time as fosolve. for the above you mentioned, Do you mean I should work on Marching squares and heuristics to speed up the calculation process like desmos? – dali May 18 '22 at 08:16
  • I would try to use a compiled version, how partial that may be. Such as with `numba.jit`. As with every dynamically typed language, there is a large overhead of type-checking and selecting appropriate variants of operations in every step. Introducing static typing at some level reduces that. Explicit loops are slow, this is the reason for `numpy`. I do not know if the `itertools` package provides a comparable acceleration. Generalist function calls have to sanitize their inputs in every call, introducing wrapper calls etc. This makes their use slow even if the initial guess is perfect. – Lutz Lehmann May 18 '22 at 08:39
  • So my guess is that a marching-squares algorithm together with some bracketing method like Dekker or Brent, implemented to make full use of JIT, could be considerably faster. Or if you have the derivatives of the function, follow the tangent vector field with an ODE method and correct using Newton at the points where you want the values. – Lutz Lehmann May 18 '22 at 08:43
  • I used to add numba.jit intro the function as decorator, but I don't how to make it works, cuz it reports error about my function. I will try to look into how make use of jit. I have the deratives of the function, in fact I use the coordinate I got to calculate the partial differential and used to integrate. So I don't know If it works with ODE and newton, Do you have any related implementation for you mentioned ODE method? cuz I am quite new to python and coding, I don't have experience about this – dali May 18 '22 at 09:43
  • For the compiled version, do you mean using Cpp language or easily trying what you meantioned packages 'numba, and itertools' – dali May 18 '22 at 09:46
  • In principle anything should work that makes the functions used strongly typed, which includes numba and cython. Some very loosely relevant examples comparing both approaches can be found in https://stackoverflow.com/questions/27346232/correctly-annotate-a-numba-function-using-jit, https://stackoverflow.com/questions/42838103/how-can-i-use-cython-well-to-solve-a-differential-equation-faster, https://stackoverflow.com/questions/26823312/numba-or-cython-acceleration-in-reaction-diffusion-algorithm – Lutz Lehmann May 18 '22 at 10:39

0 Answers0