How to compute system of multivariable equations in Python using Newton-Raphson method?

Question

I need your help with programming. Here is my system of equations:

f_1(x1, x2) = 89.3885624 + 169.6377442*x1 + 169.439564*x2 
+ 65.07923*((3*x1*(-a1-sqrt(a1^2 - 4*a0*a2)/(2*a2) +0.55071844)/3*x2) 
+ 162.698313*((-a1-sqrt(a1^2 - 4*a0*a2))/2*a2) 
+ 174.39264*x1*((3*x1*(-a1-sqrt(a1^2 - 4*a0*a2))/(2*a2) +0.55071844)/3*x2) 
+ 174.39264*x2*((-a1-sqrt(a1^2 - 4*a0*a2))/2*a2) 
+ 72.077218*x1*((-a1-sqrt(a1^2 - 4*a0*a2))/2*a2) 
+ 189.511738*x2*((3*x1*(-a1-sqrt(a1^2 - 4*a0*a2)/(2*a2) +0.55071844)/3*x2)

and

f_2(x1,x2) = 78.183644 + 26.71298*x1 + 66.782413*x2 
- 169.637744*((3*x1*(-a1-sqrt(a1^2 - 4*a0*a2)/(2*a2) +0.55071844)/3*x2) 
- 169.439564*((-a1-sqrt(a1^2 - 4*a0*a2))/2*a2) 
- 174.39264*((3*x1*(-a1-sqrt(a1^2 - 4*a0*a2)/(2*a2) +0.55071844)/3*x2)^2 
- 261.5889567*((3*x1*(-a1-sqrt(a1^2 - 4*a0*a2)/(2*a2) +0.55071844)/3*x2)
*((-a1-sqrt(a1^2 - 4*a0*a2))/2*a2) 
- 174.39264*((-a1-sqrt(a1^2 - 4*a0*a2))/2*a2)^2 
- 54.306279*x1*((-a1-sqrt(a1^2 - 4*a0*a2))/2*a2) 
+ 54.306279*x2*((3*x1*(-a1-sqrt(a1^2 - 4*a0*a2)/(2*a2) +0.55071844)/3*x2)

a0, a1 and a2 are equations with variables x1 and x2.

These two equation needs to be equal to 0 or very close to 0. I want to use Newton-Raphson method, but I do not know how. On the internet I find a lot of examples but they use easier system of equation as is mine.

Sorry for my bad English and thank you for your help!

Answer 1

You could use either scipy or mpmath, check out this answer Find roots of a system of equations to an arbitrary decimal precision

If you wanted to work it out on your own then check out this link: http://fourier.eng.hmc.edu/e176/lectures/NM/node21.html

You'll need to compute the partial derivatives of your functions to be able to compute the Jacobian matrix and then need to invert it. Your functions are pretty gnarly so you might consider using wolframscript which is a command line version of Mathematica. It'll make it a lot easier to compute the derivatives.

For example using wolframscript (with slight changes like sqrt() -> Sqrt[] and having the + at the end of the line to do multiline functions) you can easily get the partial derivatives that you need.

In[8]:= f1[x1_, x2_] := 89.3885624 + 169.6377442*x1 + 169.439564*x2 +
65.07923*((3*x1*(-a1-Sqrt[a1^2 - 4*a0*a2])/(2*a2) +0.55071844)/3*x2) +
162.698313*((-a1-Sqrt[a1^2 - 4*a0*a2])/2*a2) +
174.39264*x1*((3*x1*(-a1-Sqrt[a1^2 - 4*a0*a2])/(2*a2) +0.55071844)/3*x2) +
174.39264*x2*((-a1-Sqrt[a1^2 - 4*a0*a2])/2*a2) +
72.077218*x1*((-a1-Sqrt[a1^2 - 4*a0*a2])/2*a2) +
189.511738*x2*((3*x1*(-a1-Sqrt[a1^2 - 4*a0*a2])/(2*a2) +0.55071844)/3*x2)

In[9]:= D[f1[x1, x2], x1]                                                                          
                                           2
Out[9]= 169.638 + 36.0386 a2 (-a1 - Sqrt[a1  - 4 a0 a2]) +

                           2                                        2
     32.5396 (-a1 - Sqrt[a1  - 4 a0 a2]) x2   87.1963 (-a1 - Sqrt[a1  - 4 a0 a2]) x1 x2
>    -------------------------------------- + ----------------------------------------- +
                       a2                                        a2

                                         2
                         3 (-a1 - Sqrt[a1  - 4 a0 a2]) x1
>    58.1309 (0.550718 + --------------------------------) x2 +
                                       2 a2

                           2               2
     94.7559 (-a1 - Sqrt[a1  - 4 a0 a2]) x2
>    ---------------------------------------
                       a2