Mercurial > repos > public > sbplib
comparison +time/SBPInTimeSecondOrderFormImplicit.m @ 707:0de70ec8bf60 feature/quantumTriangles
merge with feature/optim
author | Ylva Rydin <ylva.rydin@telia.com> |
---|---|
date | Fri, 10 Nov 2017 14:22:56 +0100 |
parents | d4b999585af1 |
children | 94bd0f3293c8 |
comparison
equal
deleted
inserted
replaced
696:7c16b5af8d98 | 707:0de70ec8bf60 |
---|---|
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) | |
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 -C, -B; | |
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 |