des.math.method.calculus.integ source code

1 module des.math.method.calculus.integ;
2 
3 import des.ts;
4 import des.stdx.traits;
5 
6 import std.math;
7 import std.meta;
8 
9 ///
10 T euler(T,E1,E2,E3)( in T x, T delegate(in T,E1) f, E2 time, E3 h )
11     if( hasBasicMathOp!T && allSatisfy!(isFloatingPoint,E1,E2,E3) )
12 {
13     return x + f( x, time ) * h;
14 }
15 
16 ///
17 T runge(T,E1,E2,E3)( in T x, T delegate(in T,E1) f, E2 time, E3 h )
18     if( hasBasicMathOp!T && allSatisfy!(isFloatingPoint,E1,E2,E3) )
19 {
20     T k1 = f( x, time ) * h;
21     T k2 = f( x + k1 * 0.5, time + h * 0.5 ) * h;
22     T k3 = f( x + k2 * 0.5, time + h * 0.5 ) * h;
23     T k4 = f( x + k3, time + h ) * h;
24     return cast(T)( x + ( k1 + k2 * 2.0 + k3 * 2.0 + k4 ) / 6.0 );
25 }
26 
27 unittest
28 {
29     double a1 = 0, a2 = 0, pa = 5;
30     double time = 0, ft = 10, step = .01;
31 
32     auto rpart( in double A, double time ) { return pa; }
33 
34     foreach( i; 0 .. 1000 )
35     {
36         a1 = euler( a1, &rpart, time+=step, step );
37         a2 = runge( a1, &rpart, time+=step, step );
38     }
39 
40     auto va = ft * pa;
41 
42     assertEqApprox( va, a1, step * pa * 2 );
43     assertEqApprox( va, a2, step * pa );
44 
45     auto rpart2( in float A, double time ) { return pa; }
46 
47     static assert( !is(typeof( euler( a1, &rpart2, 0, 0 ) )));
48 }
49 
50 unittest
51 {
52     import des.math.util.mathstruct;
53 
54     static struct Pos
55     {
56         double x=0, y=0;
57         mixin( BasicMathOp!"x y" );
58     }
59 
60     static struct Point
61     {
62         Pos pos, vel;
63         mixin( BasicMathOp!"pos vel" );
64     }
65 
66     Pos acc( in Pos p )
67     { return Pos( -(p.x * abs(p.x)), -(p.y * abs(p.y)) ); }
68 
69     Point rpart( in Point p, double time )
70     { return Point( p.vel, acc(p.pos) ); }
71 
72     auto state1 = Point( Pos(50,10), Pos(5,15) );
73     auto state2 = Point( Pos(50,10), Pos(5,15) );
74 
75     double t = 0, ft = 10, dt = 0.01;
76 
77     foreach( i; 0 .. 1000 )
78     {
79         state1 = euler( state1, &rpart, t+=dt, dt );
80         state2 = runge( state2, &rpart, t+=dt, dt );
81     }
82 }
83 
84 ///
85 unittest
86 {
87     import des.math.linear.vector;
88     import des.math.util.mathstruct;
89 
90     static struct Point
91     {
92         vec3 pos, vel;
93         mixin( BasicMathOp!"pos vel" );
94     }
95 
96     Point rpart( in Point p, double time )
97     { return Point( p.vel, vec3(0,0,0) ); }
98 
99     auto v = Point( vec3(10,3,1), vec3(5,4,3) );
100 
101     double time = 0, ft = 10, step = .01;
102     foreach( i; 0 .. 1000 )
103         v = runge( v, &rpart, time+=step, step );
104 
105     assert( eq_approx( v.pos, vec3(60,43,31),  1e-3 ) );
106 }