使用scipy来解非线性方程

我们有两个正态分布：N(2,4)和N(3,5)，现在想求出靠近两者均值的等概率点，如图所示：



其中Difference=N(2,4)-N(3,5)，即为我们下列代码中的函数f

#! /usr/bin/python3
#-*-coding:utf-8-*-

from scipy.optimize import fsolve
import numpy as np
import math

#定义函数（方程f(x|arg)=0）
def f(x,*arg): #len(arg)=4, arg=(miu_1,miu_2,sigma_1,sigma_2)
f1=math.exp(-(x-arg[0])**2/2/arg[2]**2)/arg[2]/math.sqrt(2*math.pi)
f2=math.exp(-(x-arg[1])**2/2/arg[3]**2)/arg[3]/math.sqrt(2*math.pi)
return(f1-f2)

x0=(2.5,2,3,4,5)
result=fsolve(f,x0=x0[0],args=x0[1:])
print(result)
#输出：[ 5.19949597]

可以看到解出的根和我们图上的标注的位置非常接近，也就是说这个根是我们所需要的的解。

注意在fsolve函数里，得注明前多少个f的参数是我们方程中要求解的变量(x0的长度)，哪些是方程中的参数(args)

对于具有多个根的方程，x0的具体数值很重要，因为初解很大程度上决定了根收敛的方向。

为了求得靠近两分布中心位置的根，初解最好选为两均值的均值，在这里我们将x0设为(2+3)/2=2.5

如果设为1，可以看到Difference在1处的导数几乎为0，因而根会收敛得非常奇葩，实际操作做了一下，输出的结果为[-152.18013068]，远远偏离分布中心位置。因为fsolve有一个默认参数xtol=1.49012e-08，只要当两次迭代结果之间的差值小于xtol就会终止迭代，所以不难理解为什么方程只有两个不大的根，却得出了一个这么奇怪的结果。

如果将x0设为0，就可以得到另一个根：[-4.75505152]，这个是较为远离分布中心的根。

因为在这里我们只求解一个变量x，所以在定义函数的时候，只需返回一个值（一个约束），如果是求解多元方程，则在定义函数的时候返回一个list即可，fsolve会自动求使得list的所有元素均为0的根