des.math.method.calculus.diff source code

1 module des.math.method.calculus.diff;
2 
3 import des.math.linear;
4 import std.traits;
5 
6 import des.util.testsuite;
7 
8 version(unittest) import std.math;
9 
10 ///
11 auto df(size_t N, size_t M, T, E=T)
12     ( Vector!(M,T) delegate( in Vector!(N,T) ) f, in Vector!(N,T) p, E step=E.epsilon*10 )
13     if( isFloatingPoint!T && isFloatingPoint!E )
14 {
15     Matrix!(M,N,T) ret;
16 
17     T dstep = 2.0 * step;
18     foreach( i; 0 .. N )
19     {
20         Vector!(N,T) p1 = p;
21         p1[i] -= step;
22         Vector!(N,T) p2 = p;
23         p2[i] += step;
24 
25         auto r1 = f(p1);
26         auto r2 = f(p2);
27 
28         ret.setCol( i, (r2-r1) / dstep );
29     }
30 
31     return ret;
32 }
33 
34 ///
35 unittest
36 {
37     auto func( in dvec2 p ) { return dvec3( p.x^^2, sqrt(p.y) * p.x, 3 ); }
38 
39     auto res = df( &func, dvec2(18,9), 1e-5 );
40     auto must = Matrix!(3,2,double)( 36, 0, 3, 3, 0, 0 );
41 
42     assert( eq_approx( res.asArray, must.asArray, 1e-5 ) );
43 }
44 
45 ///
46 auto df_scalar(T,K,E=T)( T delegate(T) f, K p, E step=E.epsilon*2 )
47     if( isFloatingPoint!T && isFloatingPoint!E && is( K : T ) )
48 {
49     Vector!(1,T) f_vec( in Vector!(1,T) p_vec )
50     { return Vector!(1,T)( f( p_vec[0] ) ); }
51     return df( &f_vec, Vector!(1,T)(p), step )[0][0];
52 }
53 
54 ///
55 unittest
56 {
57     auto pow2( double x ){ return x^^2; }
58     auto res1 = df_scalar( &pow2, 1 );
59     auto res2 = df_scalar( &pow2, 3 );
60     auto res3 = df_scalar( &pow2, -2 );
61     assert( abs(res1 - 2.0) < 2e-6 );
62     assert( abs(res2 - 6.0) < 2e-6 );
63     assert( abs(res3 + 4.0) < 2e-6 );
64 }