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.ts;
7 
8 ///
9 auto df(size_t N, size_t M, T, E=T)
10     ( Vector!(M,T) delegate( in Vector!(N,T) ) f, in Vector!(N,T) p, E step=E.epsilon*10 )
11     if( isFloatingPoint!T && isFloatingPoint!E )
12 in { assert( f !is null ); } body
13 {
14     Matrix!(M,N,T) ret;
15     Vector!(N,T) p1, p2;
16     T dstep = 2.0 * step;
17 
18     foreach( i; 0 .. N )
19     {
20         p1 = p2 = p;
21         p1[i] -= step;
22         p2[i] += step;
23         ret.setCol( i, ( f(p2) - f(p1) ) / dstep );
24     }
25 
26     return ret;
27 }
28 
29 ///
30 unittest
31 {
32     import std.math : sqrt;
33     auto func( in dvec2 p ) { return dvec3( p.x^^2, sqrt(p.y) * p.x, 3 ); }
34 
35     auto res = df( &func, dvec2(18,9), 1e-5 );
36     auto must = Matrix!(3,2,double)( 36, 0, 3, 3, 0, 0 );
37 
38     assertEqApprox( res.asArray, must.asArray, 1e-5 );
39 }
40 
41 ///
42 auto df(T,K,E=T)( T delegate(T) f, K p, E step=E.epsilon*2 )
43     if( isFloatingPoint!T && isFloatingPoint!E && is( K : T ) )
44 {
45     alias V1 = Vector!(1,T);
46     V1 f_vec( in V1 p_vec ) { return V1( f( p_vec[0] ) ); }
47     return df( &f_vec, V1(p), step )[0][0];
48 }
49 
50 ///
51 unittest
52 {
53     auto pow2( double x ){ return x^^2; }
54     assertEqApprox( df( &pow2,  1 ), 2.0, 2e-6 );
55     assertEqApprox( df( &pow2,  3 ), 6.0, 2e-6 );
56     assertEqApprox( df( &pow2, -2 ), -4.0, 2e-6 );
57 }