【解方程】scipy.optimize.solve



2017年06月06日    Author:Guofei

文章归类: 0x55_数值计算    文章编号: 5510

版权声明:本文作者是郭飞。转载随意,但需要标明原文链接,并通知本人
原文链接:https://www.guofei.site/2017/06/06/scipyfsolve.html


fsolve求方程的解

如果要求解方程:
\(\left \{ \begin{array}{ccc} f1(u1,u2,u3)=0\\ f2(u1,u2,u3)=0\\ f3(u1,u2,u3)=0 \end{array}\right.\)

那么func这么定义:

def func(x):
  u1,u2,u3=x
  return [f1(u1,u2,u3),f2(u1,u2,u3),f3(u1,u2,u3)]

代码案例1:

import scipy.optimize as opt
import numpy as np
def f(x):
    x0,x1,x2=x
    return np.array([5*x1+3,4*x0*x0-2*np.sin(x1*x2),x1*x2-1.5])

result=opt.fsolve(f,[1,1,1])
print(result)
print(f(result))

代码案例2:

如果给了Jacobian矩阵,那么迭代速度更快
Jacobian矩阵的定义是: \(\left [ \begin{array}{ccc} \frac{\partial f1}{\partial u1}&\frac{\partial f1}{\partial u2}&\frac{\partial f1}{\partial u2}\\ \\ \frac{\partial f2}{\partial u1}&\frac{\partial f2}{\partial u2}&\frac{\partial f2}{\partial u2}\\ \\ \frac{\partial f3}{\partial u1}&\frac{\partial f3}{\partial u2}&\frac{\partial f3}{\partial u2} \end{array} \right ]\)

import scipy.optimize as opt
import numpy as np
def obj_func(x):
    x0,x1,x2=x
    return [5*x1+3,4*x0*x0-2*np.sin(x1*x2),x1*x2-1.5]

def jacobian(x):
    x0, x1, x2 = x
    return [
        [0,5,0],
        [8*x0,-2*x2*np.cos(x1*x2),-2*x1*np.cos(x1*x2)],
        [0,x2,x1]
    ]
result=opt.fsolve(obj_func,[1,1,1],fprime=jacobian)
print(result)
print(obj_func(result))

您的支持将鼓励我继续创作!