Mercurial > repos > public > sbplib

comparison +time/ExplicitRungeKuttaSecondOrderDiscreteData.m @ 858:335b8b1366d4 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add RK for second order and discrete data.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 02 Oct 2018 13:43:33 -0700
parents
children
comparison
equal deleted inserted replaced
857:3c916a00033f 858:335b8b1366d4
1 classdef ExplicitRungeKuttaSecondOrderDiscreteData < time.Timestepper
2 properties
3 k
4 t
5 w
6 m
7 D
8 E
9 M
10 C_cont % Continuous part (function handle) of forcing on first order form.
11 C_discr% Discrete part (matrix) of forcing on first order form.
12 n
13 order
14 tsImplementation % Time stepper object, RK first order form,
15 % which this wraps around.
16 end
17
18
19 methods
20 % Solves u_tt = Du + Eu_t + S by
21 % Rewriting on first order form:
22 % w_t = M*w + C(t)
23 % where
24 % M = [
25 % 0, I;
26 % D, E;
27 % ]
28 % and
29 % C(t) = [
30 % 0;
31 % S(t)
32 % ]
33 % D, E, should be matrices (or empty for zero)
34 % They can also be omitted by setting them equal to the empty matrix.
35 % S = S_cont + S_discr, where S_cont is a function handle
36 % S_discr a matrix of data vectors, one column per stage.
37 function obj = ExplicitRungeKuttaSecondOrderDiscreteData(D, E, S_cont, S_discr, k, t0, v0, v0t, order)
38 default_arg('order', 4);
39 default_arg('S_cont', []);
40 default_arg('S_discr', []);
41 obj.D = D;
42 obj.E = E;
43 obj.m = length(v0);
44 obj.n = 0;
45
46 default_arg('D', sparse(obj.m, obj.m) );
47 default_arg('E', sparse(obj.m, obj.m) );
48
49 obj.k = k;
50 obj.t = t0;
51 obj.w = [v0; v0t];
52
53 I = speye(obj.m);
54 O = sparse(obj.m,obj.m);
55
56 obj.M = [
57 O, I;
58 D, E;
59 ];
60
61 % Build C_cont
62 if ~isempty(S_cont)
63 obj.C_cont = @(t)[
64 sparse(obj.m,1);
65 S_cont(t)
66 ];
67 else
68 obj.C_cont = [];
69 end
70
71 % Build C_discr
72 if ~isempty(S_discr)
73 [~, nt] = size(S_discr);
74 obj.C_discr = [sparse(obj.m, nt);
75 S_discr
76 ];
77 else
78 obj.C_discr = [];
79 end
80 obj.tsImplementation = time.ExplicitRungeKuttaDiscreteData(obj.M, obj.C_cont, obj.C_discr,...
81 k, obj.t, obj.w, order);
82 end
83
84 function [v,t,U,T,K] = getV(obj)
85 [w,t,U,T,K] = obj.tsImplementation.getV();
86
87 v = w(1:end/2);
88 U = U(1:end/2, :); % Stage approximations in previous time step.
89 K = K(1:end/2, :); % Stage rates in previous time step.
90 % T: Stage times in previous time step.
91 end
92
93 function [vt,t,U,T,K] = getVt(obj)
94 [w,t,U,T,K] = obj.tsImplementation.getV();
95
96 vt = w(end/2+1:end);
97 U = U(end/2+1:end, :); % Stage approximations in previous time step.
98 K = K(end/2+1:end, :); % Stage rates in previous time step.
99 % T: Stage times in previous time step.
100 end
101
102 function [a,b,c,s] = getTableau(obj)
103 [a,b,c,s] = obj.tsImplementation.getTableau();
104 end
105
106 % Returns quadrature weights for stages in one time step
107 function quadWeights = getTimeStepQuadrature(obj)
108 [~, b] = obj.getTableau();
109 quadWeights = obj.k*b;
110 end
111
112 % Use RK for first order form to step
113 function obj = step(obj)
114 obj.tsImplementation.step();
115 [v, t] = obj.tsImplementation.getV();
116 obj.w = v;
117 obj.t = t;
118 obj.n = obj.n + 1;
119 end
120 end
121
122 methods (Static)
123 function k = getTimeStep(lambda, order)
124 default_arg('order', 4);
125 k = obj.tsImplementation.getTimeStep(lambda, order);
126 end
127 end
128
129 end