Как вычислить систему уравнений с несколькими переменными в Python, используя метод Ньютона-Рафсона?

Мне нужна ваша помощь в программировании. Вот моя система уравнений:

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)

а также

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 и a2 — уравнения с переменными x1 и x2.

Эти два уравнения должны быть равны 0 или очень близки к 0. Я хочу использовать метод Ньютона-Рафсона, но не знаю как. В Интернете я нахожу много примеров, но они используют более простую систему уравнений, как моя.

Извините за мой плохой английский и спасибо за вашу помощь!


python function math equation newtons-method
person julio77    schedule 09.09.2020    source источник
comment
Посетите эту страницу fourier.eng.hmc.edu/e176/lectures/ NM/node21.html в нем есть несколько примеров системы уравнений, похожей на то, что вы описываете.   -  person g23    schedule 10.09.2020

Ответы (1)


arrow_upward
1
arrow_downward

Вы можете использовать либо scipy, либо mpmath, посмотрите этот ответ Найти корни системы уравнений с произвольной десятичной точностью

Если вы хотите решить это самостоятельно, перейдите по этой ссылке: http://fourier.eng.hmc.edu/e176/lectures/NM/node21.html

Вам нужно будет вычислить частные производные ваших функций, чтобы иметь возможность вычислить матрицу Якоби, а затем инвертировать ее. Ваши функции довольно корявые, поэтому вы можете рассмотреть возможность использования wolframscript, который представляет собой версию Mathematica для командной строки. Это значительно облегчит вычисление производных.

Например, используя wolframscript (с небольшими изменениями, такими как sqrt() -> Sqrt[] и наличие + в конце строки для выполнения многострочных функций), вы можете легко получить нужные вам частные производные.

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
person g23    schedule 11.09.2020