Mercurial > repos > public > sbplib

comparison +time/SBPInTimeSecondOrderFormImplicit.m @ 886:8894e9c49e40 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default for latest changes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 15 Nov 2018 16:36:21 -0800
parents +time/SBPInTimeSecondOrderForm.m@b5e5b195da1e +time/SBPInTimeSecondOrderForm.m@66eb4a2bbb72
children f5e14e5986b5
comparison
equal deleted inserted replaced
816:b5e5b195da1e 886:8894e9c49e40
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, do_scaling, TYPE, order, blockSize)
18 default_arg('f', []);
19 default_arg('TYPE', []);
20 default_arg('order', []);
21 default_arg('blockSize',[]);
22 default_arg('do_scaling', false);
23
24 m = length(v0);
25
26 default_arg('A', speye(m, m));
27 default_arg('B', sparse(m, m));
28 default_arg('C', sparse(m, m));
29
30 I = speye(m);
31 O = sparse(m,m);
32
33 % Rewrite to
34 % AA*w_t = BB*w + ff(t);
35
36 obj.AA = [
37 I, O;
38 O, A;
39 ];
40 obj.BB = [
41 O, I;
42 -C, -B;
43 ];
44
45 if ~isempty(f)
46 obj.ff = @(t)[
47 sparse(m,1);
48 f(t);
49 ];
50 else
51 obj.ff = @(t) sparse(2*m,1);
52 end
53
54 w0 = [v0; v0t];
55
56 obj.k = k;
57 obj.t = t0;
58 obj.n = 0;
59
60 if do_scaling
61 scaling = [ones(m,1); sqrt(diag(C))];
62 obj.firstOrderTimeStepper = time.SBPInTimeScaled(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, scaling, TYPE, order, blockSize);
63 else
64 obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
65 end
66 end
67
68 function [v,t] = getV(obj)
69 w = obj.firstOrderTimeStepper.getV();
70 v = w(1:end/2);
71 t = obj.t;
72 end
73
74 function [vt,t] = getVt(obj)
75 w = obj.firstOrderTimeStepper.getV();
76 vt = w(end/2+1:end);
77 t = obj.t;
78 end
79
80 function state = getState(obj)
81 [v, t] = obj.getV();
82 [vt] = obj.getVt();
83 state = struct('v', v, 'vt', vt, 't', t, 'k', obj.k);
84 end
85
86 function obj = step(obj)
87 obj.firstOrderTimeStepper.step();
88 obj.t = obj.firstOrderTimeStepper.t;
89 obj.n = obj.firstOrderTimeStepper.n;
90 end
91 end
92 end