Mercurial > repos > public > sbplib

annotate +time/ExplicitRungeKuttaDiscreteData.m @ 753:44c46bd6913a feature/d1_staggered

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