python样条插值的实现（四）

上一次我们完成了对第二类边界条件的样条插值的函数编写，现在的问题是：如何求函数上的最大曲率

以我非常有限的python知识，想到了三种办法：

1.先将曲率的表达式给出，再给足x的采样点，得到x对应的曲率的值列表，利用列表的寻找最大值函数来找曲率最大值

当然这个方法弊端比较明显，首先是精度问题，精度肯定不是很高，并且计算量也较大，这里就不给出代码了（因为被我删了）

当然也可以利用取样比较再迭代的方法求得高精度，不过这是后话。

不过max函数有个小技巧，例如一个嵌套列表：[[1,50],[2,40],[5,100]]，直接用max函数肯定是不能对其进行排序的，但我们可以：

res = max(cur_all, key=lambda x: x[1])

利用lambda函数就可以得到列表中第1个（列表从第0个开始）的最大值所在的列表，即返回值应当是

res=[5,100]

 

2.将曲率表达式得出后，利用函数求极值的方法求得最大值

首先我们知道，在样条插值函数里，曲率是连续的，导数分段连续，所以我们可以根据每段的导数正负性找到每段的最大值，再用列表求最大值函数找到全局最大值。

当然，这个方法的优缺点也很明显，优点是给出的最优值精确度较高，计算量相对也不是太大。缺点是实现起来较为复杂。

3.利用scipy库中的求局部最小值函数

官方文档：

scipy.optimize.minimize_scalar(fun, bracket=None, bounds=None, args=(), method='brent', tol=None, options=None)

举个栗子：

from scipy.optimize import minimize_scalar

def f(x):
     return (x - 2) * x * (x + 2)**2

res = minimize_scalar(f)
res.x
# 1.28077640403

# 利用bounded方法，可以找指定范围内的最小值
res = minimize_scalar(f, bounds=(-3, -1), method='bounded')
res.x
# -2.0000002026

最后我采用了这个方法来求最大值。什么？你说只能求最小值？函数取个负号就是最大值了嘛！

上代码：

def curvature(f,xmin,xmax):
    # 格式为sympy函数表达式
    t=s.symbols('t')
    df=s.diff(f,t)
    ddf=s.diff(df,t)
    k=-abs(ddf)/(1+df**2)**1.5
    kk = s.lambdify(t, k)
    m = opt.minimize_scalar(kk, bounds=(xmin,xmax), method='bounded') # 取最小值
    # print(m)
    return [m.x,-m.fun]

传入函数是sympy的symbols符号函数，返回值时曲率的最大值以及最大值的横坐标。

接下来就是找到最小的两个首尾一阶导数值使得样条插值的曲率最大值最小。

由于直接通过两个一阶导数求曲率最大值的理论推导有点复杂，所以我在这里采用的是梯度下降法求，代码如下：

step=1e-6;eta=5
diff_start=-0;diff_end=9
last_diff=0;res_y1=[0,0]
diff_x=[];diff_y=[]
while True:
    last_diff=res_y1[1]
    fx, res_x1 = cur_all_group(x, y, diff_start - step, diff_end)
    fx, res_x2 = cur_all_group(x, y, diff_start + step, diff_end)
    fx, res_y1 = cur_all_group(x, y, diff_start, diff_end - step)
    fx, res_y2 = cur_all_group(x, y, diff_start, diff_end + step)
    dx=(res_x2[1]-res_x1[1])/(2*step)
    dy=(res_y2[1]-res_y1[1])/(2*step)
    diff_start=diff_start-dx*eta
    diff_end=diff_end-dy*eta
    print(dx,'dx',dy,'dy',len(diff_x)+1)
    print(res_y1[1],'cur')
    print(diff_start,diff_end,'start,end')
    diff_x.append(diff_start)
    diff_y.append(diff_end)
    if(abs(res_y1[1]-last_diff)<1e-5):
        break
my_xi=[]
my_yi=[]
t=s.symbols('t')
for i in range(len(fx)):
    xi_g=np.arange(x[i],x[i+1],0.1)
    for j in range(len(xi_g)):
        my_yi.append(fx[i].evalf(subs={t:xi_g[j]}))
        my_xi.append(xi_g[j])
plt.plot(my_yi,my_xi,'b',label = 'spline best')

在两个方向通过算出导数值，通过学习率来调整下降速度，逐步迭代找到局部最优解。

初始值为0,9时（初始值可以通过前两个点连线的斜率以及后两个点连线的斜率确定），迭代了46次得到了局部最优解：

看图中的效果确实比前两个更好一些。

全部的代码文件我打包放在资源库里，有需要自行下载。