Mercurial > repos > public > sbplib

comparison +time/SBPInTimeSecondOrderFormImplicit.m @ 460:e0caae9ef6ed feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add SBPinTime for linear DAE formulations (BUGS!)
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 19 Jun 2017 16:50:13 +0200
parents +time/SBPInTimeSecondOrderForm.m@c92d2f8319c2
children 45a3924140f4
comparison
equal deleted inserted replaced
459:1147db8a2ffa 460:e0caae9ef6ed
1 classdef SBPInTimeSecondOrderFormImplicit < time.Timestepper
2 properties
3 A, B, C, f
4 AA, BB, ff
5
6 n
7 t
8 k
9
10 firstOrderTimeStepper
11 end
12
13 methods
14 % Solves A*u_tt + B*u_t + C*u + f(t) = 0
15 % A, B can either both be constants or both be function handles,
16 % They can also be omitted by setting them equal to the empty matrix.
17 function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, TYPE, order, blockSize)
18 default_arg('f', []);
19 default_arg('TYPE', []);
20 default_arg('order', []);
21 default_arg('blockSize',[]);
22
23 m = length(v0);
24
25 default_arg('A', sparse(m, m));
26 default_arg('B', sparse(m, m));
27 default_arg('C', sparse(m, m));
28
29 I = speye(m);
30 O = sparse(m,m);
31
32 % Rewrite to
33 % AA*w_t = BB*w + ff(t);
34
35 obj.AA = [
36 I, O;
37 O, A;
38 ];
39 obj.BB = [
40 O, I;
41 -B, -C;
42 ];
43
44 if ~isempty(f)
45 obj.ff = @(t)[
46 sparse(m,1);
47 -f(t);
48 ];
49 else
50 obj.ff = @(t) sparse(2*m,1);
51 end
52
53 w0 = [v0; v0t];
54
55 obj.k = k;
56 obj.t = t0;
57 obj.n = 0;
58
59 obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
60 end
61
62 function [v,t] = getV(obj)
63 w = obj.firstOrderTimeStepper.getV();
64 v = w(1:end/2);
65 t = obj.t;
66 end
67
68 function [vt,t] = getVt(obj)
69 w = obj.firstOrderTimeStepper.getV();
70 vt = w(end/2+1:end);
71 t = obj.t;
72 end
73
74 function obj = step(obj)
75 obj.firstOrderTimeStepper.step();
76 obj.t = obj.firstOrderTimeStepper.t;
77 obj.n = obj.firstOrderTimeStepper.n;
78 end
79 end
80 end