0

I am trying to solve a complex non-linear system. The problem is that there are many roots, and some of them are complex roots. From these roots I need to select only the one that the real part is between [0,1] and does not have complex part (ex: 0.23+0i). ex:

root1: 1.02 + 2i

root2: 0.23 + 1.23i

root3: 0.23 + 0i

...

Here is my system: tau1 and tau2 are the variables that I need to find. The equations are t1 and t2 that are dependent on tau1 and tau2

x0=0 # initial position
xf=30 # final position
x1= 10;
x2 = 20;
tf=20 # final time

tau_wp=[]
def f(tau_wp):
    tau1, tau2 = tau_wp
    a=(1-tau1)**5*(-10*tau2**3 +15*tau2**4 -6*tau2**5) + (1-tau1)**4*(20*tau2**3 -35*tau2**4 + 15*tau2**5) + (1-tau1)**3*(-10*tau2**3 +20*tau2**4 -10*tau2**5) + (tau2-tau1)**5
    b=(1-tau2)**5*(-10*tau1**3 +15*tau1**4 -6*tau1**5) + (1-tau2)**4*(20*tau1**3 -35*tau1**4 + 15*tau1**5) + (1-tau2)**3*(-10*tau1**3 +20*tau1**4 -10*tau1**5)

    den=a*b -36*tau2**5*(1-tau2)**5*tau1**5*(1-tau1)**5

    p2=(-6*tau1**5*(1-tau1)**5*(xf-x0)*(10*tau2**3-15*tau2**4+6*tau2**5) \
                                         -(xf-x0)*(10*tau1**3-15*tau1**4+6*tau1**5)*a \
                                         + (x1-x0)*a + (x2-x0)*(6*tau1**5*(1-tau1**5)))

    p1=(a*b -36*tau2**5*(1-tau2)**5*tau1**5*(1-tau1)**5)*(-(xf-x0)*(10*tau1**3-15*tau1**4+6*tau1**5) +(x1-x0)) \
                        + b*( (xf-x0)*(10*tau2**3 -15*tau2**4 +6*tau2**5)*(6*tau1**5*(1-tau1)**5) \
                        + (xf-x0)*(10*tau1**3 -15*tau1**4 +6*tau1**5)*a \
                        - (x1-x0)*a - (x2-x0)*(6*tau1**5*(1-tau1)**5))


    u0=(xf-x0)*(30*tau1**2-60*tau1**3 +30*tau1**4)+p1*tf**5/120*(60*tau1**9-270*tau1**8+480*tau1**7-420*tau1**6+180*tau1**5-30*tau1**4) \
           + p2*tf**5/120 * ((1-tau2)**5*(-30*tau1**2 +60*tau1**3 -30*tau1**4) + (1-tau2)**4*(60*tau1**2 - 140*tau1**3 +75*tau1**4) + \
            (1-tau2)**3*(-30*tau1**2 +80*tau1**3 - 50*tau1**4))

    u1=(xf-x0)*(30*tau2**2 - 60*tau2**3 + 30*tau2**4)+p1*tf**5/120*((1-tau1)**5*(-30*tau2**2 +60*tau1**3 -30*tau1**4) + \
        (1-tau1)**4*(60*tau2**2 - 140*tau2**3 +75*tau2**4) + (1-tau1)**3*(-30*tau2**2 +80*tau2**3 - 50*tau2**4) \
        + 5*tau2**4 -20*tau2**3*tau1 +30*tau2**2*tau1**2 -20*tau2*tau1**3 +5*tau1**4) \
        + p2*tf**5/120*(60*tau2**9-270*tau2**8+480*tau2**7-420*tau2**6+180*tau2**5-30*tau2**4)

    ## system of nonlinear equations dependent on tau1 and tau2
    t1=u0*p1 ### equation 1
    t2=u1*p2 ### equation 2

    return [t1,t2]

I tried to use fsolve, but with fsolve I couldn't get the complex part.

Is there any way to do this in python?

Thank you so much for your help!

João Pereira
  • 25
  • 4
  • You need to be more specific. For an arbitrary nonlinear system this is impossible: there might be infinitely many roots and they might not be expressible in any closed form. It looks like you want to find the real roots of an univariate equation. Is it a polynomial? If so sympy's `real_roots` function might be what you want. Otherwise you'll have to specify something about the kind of equation that it is – Oscar Benjamin Apr 13 '20 at 15:14
  • @OscarBenjamin , thank you for your answer. I just updated the problem. I'd be grateful to have your help once again. – João Pereira Apr 13 '20 at 15:37
  • 1
    This is a polynomial system and is extremely ill-conditioned so you should not approach this using a numerical library. Take a look at sympy – Oscar Benjamin Apr 13 '20 at 17:35

1 Answers1

0

It's a little confusing because you say that you "couldn't get the complex part" but in the question you say that you are looking for solutions where the imaginary part is 0 and the magnitude of the real part is between 0 and 1. If this is correct, then nsolve can solve this pair of equations if you give a good enough initial guess:

>>> from sympy import symbols
>>> v = symbols('tau1:3')
>>> nsolve(f(v), (tau1, tau2), (.5,.4))
Matrix([
[0.495387590772031],
[ 0.49736468918969]])

You can get a rough idea of where to look for roots by looking at the values of t1 and t2 for different values of tau1 and tau2. Since they should both be zero I look at the log of the sum of squares -- the smaller the better:

>>> Matrix(10,10,lambda i,j:
log(sqrt(sum([k.subs(tau1,i/10).subs(tau2,j/10)**2 for k in (t1,t2)]))).round())
Matrix([
[zoo, zoo, zoo, zoo, zoo, zoo, zoo, zoo, zoo, zoo],
[-17,  -5,  -6,  -3,  -2,  -2,  -3,  -4,  -5,  -7],
[ -7,  -2,  -1,   0,   0,   0,   0,   0,  -1,  -2],
[ -1,   0,   3,   4,   3,   3,   4,   4,   3,   1],
[  3,   4,   6,   7,   6,   4,   7,   7,   6,   4],
[  5,   6,   8,   9,   9,   5,   9,   9,   8,   6],
[  8,   8,  10,  11,  11,   6,  11,  11,  10,   8],
[  9,   9,  11,  12,  12,   7,  12,  12,  11,   9],
[  9,  10,  12,  13,  13,   7,  13,  13,  12,  10],
[  8,  11,  13,  13,  13,   7,  13,  13,  12,  10]])

(The zoo values correspond to the trivial solution when tau1 is zero.

smichr
  • 16,948
  • 2
  • 27
  • 34