Mercurial > repos > public > sbplib

comparison +time/ExplicitRungeKuttaDiscreteData.m @ 690:9c3c180c7ed3 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement RK4 for discrete data.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 06 Mar 2018 11:59:43 -0800
parents
children 521076a6ffa5
comparison
equal deleted inserted replaced
682:ec0ac76006e2 690:9c3c180c7ed3
1 classdef ExplicitRungeKuttaDiscreteData < time.Timestepper
2 properties
3 D
4 S % Function handle for time-dependent data
5 data % Matrix of data vectors, one column per stage
6 F
7 k
8 t
9 v
10 m
11 n
12 order
13 a, b, c, s % Butcher tableau
14 K % Stage rates
15 U % Stage approximations
16 T % Stage times
17 end
18
19
20 methods
21 function obj = ExplicitRungeKuttaDiscreteData(D, S, data, k, t0, v0, order)
22 default_arg('order', 4);
23 default_arg('S', []);
24 default_arg('data', []);
25
26 obj.D = D;
27 obj.S = S;
28 obj.k = k;
29 obj.t = t0;
30 obj.v = v0;
31 obj.m = length(v0);
32 obj.n = 0;
33 obj.order = order;
34 obj.data = data;
35
36 switch order
37 case 4
38 [obj.a, obj.b, obj.c, obj.s] = time.rkparameters.rk4();
39 otherwise
40 error('That RK method is not available');
41 end
42
43 obj.K = sparse(obj.m, obj.s);
44 obj.U = sparse(obj.m, obj.s);
45
46 end
47
48 function [v,t,U,T] = getV(obj)
49 v = obj.v;
50 t = obj.t;
51 U = obj.U; % Stage approximations in previous time step.
52 T = obj.T; % Stage times in previous time step.
53 end
54
55 function [a,b,c,s] = getTableau(obj)
56 a = obj.a;
57 b = obj.b;
58 c = obj.c;
59 s = obj.s;
60 end
61
62 function obj = step(obj)
63 v = obj.v;
64 a = obj.a;
65 b = obj.b;
66 c = obj.c;
67 s = obj.s;
68 S = obj.S;
69 dt = obj.k;
70 K = obj.K;
71 U = obj.U;
72 D = obj.D;
73 data = obj.data;
74
75 for i = 1:s
76 U(:,i) = v;
77 for j = 1:i-1
78 U(:,i) = U(:,i) + dt*a(i,j)*K(:,j);
79 end
80
81 K(:,i) = D*U(:,i);
82 obj.T(i) = obj.t + c(i)*dt;
83
84 % Data from continuos function and discrete time-points.
85 if ~isempty(S)
86 K(:,i) = K(:,i) + S(obj.T(i));
87 end
88 if ~isempty(data)
89 K(:,i) = K(:,i) + data(:,obj.n*s + i);
90 end
91
92 end
93
94 obj.v = v + dt*K*b;
95 obj.t = obj.t + dt;
96 obj.n = obj.n + 1;
97 obj.U = U;
98 end
99 end
100
101
102 methods (Static)
103 function k = getTimeStep(lambda)
104
105 switch obj.order
106 case 4
107 k = rk4.get_rk4_time_step(lambda);
108 otherwise
109 error('Time-step function not available for this order');
110 end
111 end
112 end
113
114 end