着手查出牛顿法不稳定的原因

今天把测试Hessian的代码加进去了，但是发现Hessian并无大的问题。

测试例子:一根圆柱向指定方向弯曲。

下面对比，在有弯曲方向参数，与无方向参数（该参数被被乘以0）的情况下，各帧迭代的函数值。

无方向参数:

帧号	1	2	3	4	5
第1次迭代	0	73777.396480	615269.331763	1566548.018724	2799725.912609
第2次迭代		150.018686	1116.876970	3499.198322	8757.705770
第3次迭代		144.756549	785.078791	1510.680152	2220.662691
第4次迭代			785.053385	1510.524488	2219.948136

有方向参数

帧号	1	2	3	4	5
第1次迭代	0	73777.396480	615457.873791	1566753.791283	2797987.408846
第2次迭代		238.365110	1093.703963	2470.166542	6643.182270
第3次迭代		234.873881	837.922487	693.717725	561.082432
第4次迭代		234.873881	837.922027	693.695513	560.826661

关键在于第二帧，两种方法的初始函数值是相同的，但是最终却收敛到不同的值。而且有方向的参数收敛到的值要大一些。猜测是收敛到局部极值点了。为了验证这一点，需要检查第二帧两种方法的梯度。

	1	2	3	4
无方向梯度Gn	1437015.459672	12004.873928	9.966298
有方向梯度Gn	1437015.459672	10015.666457	1.390628	0.000000
无方向梯度Gp	5596204.567079	2756.358209	0.697293
有方向梯度Gp	5596204.567079	10439.419591	464.339724	0.271374

两个梯度都收敛，但估计收敛的方向不同。

若改变导数有限差商的步长，不影响收敛趋势，只是具体每步迭代的梯度值会稍微变化。

在计算导数的过程中，q,v,a,p的值都没有被错误修改。因此迭代的过程是确实收敛的。

写一个测试导数的命令。

用无方向的方法模拟一次，测试第二帧最终梯度:

funVal:144.757682 gradient: |gn| = 5.132273, |gp| = 1.245470

测试梯度是否为函数线性近似:
// step = 1.000000e+002 funVal = 1.525598e+014 approxVal = -2.973997e+002 error = 1.525598e+014 error/dx^2 = 1.525598e+010 //
// step = 1.000000e+001 funVal = 1.498637e+012 approxVal = 1.005419e+002 error = 1.498637e+012 error/dx^2 = 1.498637e+010 //
// step = 1.000000e+000 funVal = 1.244686e+010 approxVal = 1.403361e+002 error = 1.244686e+010 error/dx^2 = 1.244686e+010 //
// step = 1.000000e-001 funVal = 6.226740e+007 approxVal = 1.443155e+002 error = 6.226725e+007 error/dx^2 = 6.226725e+009 //
// step = 1.000000e-002 funVal = 6.349599e+005 approxVal = 1.447135e+002 error = 6.348152e+005 error/dx^2 = 6.348152e+009 //
// step = 1.000000e-003 funVal = 6.504549e+003 approxVal = 1.447533e+002 error = 6.359796e+003 error/dx^2 = 6.359796e+009 //
// step = 1.000000e-004 funVal = 2.083666e+002 approxVal = 1.447572e+002 error = 6.360932e+001 error/dx^2 = 6.360932e+009 //

此时把系统改为有方向，测试此时梯度:

funVal2:144.757682 gradient2: |gn| = 5.132273, |gp| = 1.285636

测试梯度是否为函数线性近似:

// step = 1.000000e+002 funVal = 1.412399e+014 approxVal = 2.651699e+002 error = 1.412399e+014 error/dx^2 = 1.412399e+010 //
// step = 1.000000e+001 funVal = 1.386940e+012 approxVal = 1.567989e+002 error = 1.386940e+012 error/dx^2 = 1.386940e+010 //
// step = 1.000000e+000 funVal = 1.000263e+010 approxVal = 1.459618e+002 error = 1.000263e+010 error/dx^2 = 1.000263e+010 //
// step = 1.000000e-001 funVal = 5.004513e+007 approxVal = 1.448781e+002 error = 5.004499e+007 error/dx^2 = 5.004499e+009 //
// step = 1.000000e-002 funVal = 4.971193e+005 approxVal = 1.447697e+002 error = 4.969745e+005 error/dx^2 = 4.969745e+009 //
// step = 1.000000e-003 funVal = 5.112932e+003 approxVal = 1.447589e+002 error = 4.968173e+003 error/dx^2 = 4.968173e+009 //
// step = 1.000000e-004 funVal = 1.944345e+002 approxVal = 1.447578e+002 error = 4.967674e+001 error/dx^2 = 4.967674e+009 //
// step = 1.000000e-005 funVal = 1.452514e+002 approxVal = 1.447577e+002 error = 4.937210e-001 error/dx^2 = 4.937210e+009 //

可见此时结果比较接近。现在怀疑有方向和无方向这两种情况，在迭代的后期对函数值仅有微小的影响，在迭代的后期，步长的前进极为缓慢，因而导致不同的结果。

进一步提高迭代精度发现，在迭代的后期，梯度值可以变得很小。这暗示函数是收敛到局部最小值了。