des.math.method.calculus.integ source code

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