- 论坛徽章:
- 1
|
本帖最后由 lost_templar 于 2013-04-12 17:11 编辑
自己写了做自动数值求导的;
一次导数正常,二次就会出错,有些莫名其妙。
一次求导正确:- double ds( double x )
- {
- return std::sin(x);
- }
复制代码 求导使用:- derivative<0, double(double)>( ds )( 0.0 );
复制代码 输出 1
二次求导则会出错:- auto const& dds = derivative<0, double(double)>(ds);
- derivative<0, double(double)>(dds)( 0.0 );
复制代码 依旧输出 1,不知为何
derivitative 实现的大致代码如下:- template<std::size_t M, typename Dummy_Type> struct derivative;
- template< std::size_t M, typename R, typename... Types >
- struct derivative< M, R(Types...) >
- {
- typedef R return_type;
- typedef std::function< R(Types...) > function_type;
- typedef typename type_at< M, Types... >::result_type result_type;
- typedef std::size_t size_type;
- function_type ff;
- template< typename Function >
- derivative( const Function& ff_ ) : ff( ff_ ) { }
- result_type operator()( Types... vts ) const
- {
- typedef typename stepper_value_type<result_type>::value_type value_type; //fix for complex
- const value_type decrease_step = 1.6180339887498948482;
- const value_type decrease_step_2 = 2.6180339887498948482;
- const value_type safe_boundary = 2.6180339887498948482;
- const value_type start_step = 1.0e-8;
- const result_type x = value_at<M, Types...>()( vts... );
- const size_type iter_depth = 64;
- result_type error = std::numeric_limits<value_type>::max();
- result_type ans = std::numeric_limits<value_type>::quiet_NaN();
- value_type step = start_step;
- auto const& lhs_f = [&step]( result_type x ) { return x - step; };
- auto const& rhs_f = [&step]( result_type x ) { return x + step; };
- oscillate_function< M, return_type, Types... > const lhs_of( ff, lhs_f );
- oscillate_function< M, return_type, Types... > const rhs_of( ff, rhs_f );
- //matrix<result_type> a( iter_depth, iter_depth );
- result_type a[iter_depth][iter_depth];
- a[0][0] = ( rhs_of( vts... ) - lhs_of( vts... ) ) / ( step + step );
- //for ( size_type i = 1; i != a.row(); ++i )
- for ( size_type i = 1; i != iter_depth; ++i )
- {
- step /= decrease_step;
- a[i][0] = ( rhs_of( vts... ) - lhs_of( vts... ) ) / ( step + step );
- result_type factor = decrease_step_2;
- for ( size_type j = 1; j <= i; ++j )
- {
- const result_type factor_1 = factor - result_type(1);
- a[i][j] = ( a[i][j-1] * factor - a[i-1][j-1] ) / factor_1;
- factor *= decrease_step_2;
- const result_type error_so_far = std::max( std::abs(a[i][j]-a[i][j-1]), std::abs(a[i][j]-a[i-1][j-1]) );
- if ( error > error_so_far )
- {
- error = error_so_far;
- ans = a[i][j];
- }
- }
- if ( std::abs( a[i][i] - a[i-1][i-1] ) >= safe_boundary * error )
- return ans;
- }
- return ans;
- }
- };
复制代码 其中 type_at 如下实现:- template< std::size_t N, typename T, typename... Types >
- struct type_at
- {
- static_assert( N < sizeof...(Types)+1, "dim size exceeds limitation!" );
- typedef typename type_at<N-1, Types...>::result_type result_type;
- };
- template<typename T, typename... Types>
- struct type_at< 0, T, Types...>
- {
- typedef T result_type;
- };
复制代码 那个 matrix 可以用个二维数组代替,应该不难理解。
|
|