Python 求多重根的方法

Python有不少求根的办法,但是都是求一重根的,不太适用。
在这个问题上用的是暴力搜索法:记 $$F(x)=f(x)-y$$
判断每个 x 对应的 $F(x)$ 的符号,如果$F(x_i)<0,\ F(x_{i+1})>0$,或者$F(x_i)>0,\ F(x_{i+1})<0$,则取该$F(x)=0$的根为$$x_0=(x_i+x_{i+1})/2$$

当然,这是一种比较简单粗略的方法,因为我们知道了单根所在的区间 $[x_i,\ x_{i+1}]$ 之后,就可以采用更精确的求根方法了,比如二分法、Newton-Raphson方法等等。

两个例子:
$$f_1(x)=(x-1)^2-4$$
$$f_2(x)=cos(x)$$
求根区间为 [-10,4],步长为 0.001
以上算法求出来的解为:
$root_1=-0.9995,\ 3.0005$
$root_2/\pi=-2.49984669,\ -1.50003534,\ -0.49990568,\ 0.49990568$
真值为:
$Root_1=-1,\ 3$
$Root_2/\pi=-2.5,\ -1.5,\ -0.5,\ 0.5$

代码:

from numpy import *

a=-10
b=4
d=0.001

x = arange(a, b, d)

f1 = lambda x:(x-1)**2-4
interp1 = diff(sign(f1(x))).nonzero()[0]
root1 = (x[interp1+1]+x[interp1])/2

f2 = cos(x)
interp2 = diff(sign(f2)).nonzero()[0]
root2 = (x[interp2+1]+x[interp2])/2/pi

sign() 函数为判断输入值的符号,这样区间内没有异号的话,diff(sign(f1(x))) 就为 0,有异号的情况diff(sign(f1(x))) 为 2 或者 -2。取出不为零的值则等于取出了根所在的区间即 $[x_i,\ x_{i+1}]$ 。

标签: scipy

赞 (13)

已有 2 条评论

  1. Jing Jing

    沙发党

    1. 恭喜!你是博客改版以来的第一个评论者!

添加新评论