Mercurial > repos > public > sbplib

comparison +time/CdiffImplicit.m @ 978:1a30dbe99c7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor CdiffImplicit to take input arguments in the right order
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 07 Jan 2019 16:15:49 +0100
parents ea8fefc326b2
children 47e86b5270ad
comparison
equal deleted inserted replaced
977:c9009d5a3101 978:1a30dbe99c7c
1 classdef CdiffImplicit < time.Timestepper 1 classdef CdiffImplicit < time.Timestepper
2 properties 2 properties
3 A, B, C, G 3 A, B, C
4 AA, BB, CC 4 AA, BB, CC
5 G
5 k 6 k
6 t 7 t
7 v, v_prev 8 v, v_prev
8 n 9 n
9 10
11 L,U,p,q 12 L,U,p,q
12 end 13 end
13 14
14 methods 15 methods
15 % Solves 16 % Solves
16 % A*u_tt + B*u + C*u_t = G(t) 17 % A*v_tt + B*v_t + C*v = G(t)
17 % u(t0) = f1 18 % v(t0) = v0
18 % u_t(t0) = f2 19 % v_t(t0) = v0t
19 % starting at time t0 with timestep k 20 % starting at time t0 with timestep
20 21 % Using
21 % TODO: Fix order of matrices 22 % A*Dp*Dm*v_n + B*D0*v_n + C*I0*v_n = G(t_n)
22 function obj = CdiffImplicit(A, B, C, G, f1, f2, k, t0) 23 function obj = CdiffImplicit(A, B, C, G, v0, v0t, k, t0)
23 default_arg('A', []); 24 m = length(v0);
24 default_arg('C', []); 25 default_arg('A', speye(m));
25 default_arg('G', []); 26 default_arg('B', sparse(m,m));
26 default_arg('f1', 0); 27 default_arg('G', @(t) sparse(m,1));
27 default_arg('f2', 0);
28 default_arg('t0', 0); 28 default_arg('t0', 0);
29
30 m = size(B,1);
31
32 if isempty(A)
33 A = speye(m);
34 end
35
36 if isempty(C)
37 C = sparse(m,m);
38 end
39
40 if isempty(G)
41 G = @(t) sparse(m,1);
42 end
43
44 if isempty(f1)
45 f1 = sparse(m,1);
46 end
47
48 if isempty(f2)
49 f2 = sparse(m,1);
50 end
51 29
52 obj.A = A; 30 obj.A = A;
53 obj.B = B; 31 obj.B = B;
54 obj.C = C; 32 obj.C = C;
55 obj.G = G; 33 obj.G = G;
56 34
57 AA = 1/k^2*A + 1/2*B + 1/(2*k)*C; 35 % Rewrite as AA*v_(n+1) + BB*v_n + CC*v_(n-1) = G(t_n)
58 BB = -2/k^2*A; 36 AA = A/k^2 + B/(2*k) + C/2;
59 CC = 1/k^2*A + 1/2*B - 1/(2*k)*C; 37 BB = -2*A/k^2;
60 % AA*v_next + BB*v + CC*v_prev == G(t_n) 38 CC = A/k^2 - B/(2*k) + C/2;
61 39
62 obj.AA = AA; 40 obj.AA = AA;
63 obj.BB = BB; 41 obj.BB = BB;
64 obj.CC = CC; 42 obj.CC = CC;
65 43
66 v_prev = f1; 44 v_prev = v0;
67 I = speye(m); 45 I = speye(m);
68 % v = (1/k^2*A)\((1/k^2*A - 1/2*B)*f1 + (1/k*I - 1/2*C)*f2 + 1/2*G(0)); 46 v = v0 + k*v0t;
69 v = f1 + k*f2;
70
71 47
72 if ~issparse(A) || ~issparse(B) || ~issparse(C) 48 if ~issparse(A) || ~issparse(B) || ~issparse(C)
73 error('LU factorization with full pivoting only works for sparse matrices.') 49 error('LU factorization with full pivoting only works for sparse matrices.')
74 end 50 end
75 51
77 53
78 obj.L = L; 54 obj.L = L;
79 obj.U = U; 55 obj.U = U;
80 obj.p = p; 56 obj.p = p;
81 obj.q = q; 57 obj.q = q;
82
83 58
84 obj.k = k; 59 obj.k = k;
85 obj.t = t0+k; 60 obj.t = t0+k;
86 obj.n = 1; 61 obj.n = 1;
87 obj.v = v; 62 obj.v = v;
92 v = obj.v; 67 v = obj.v;
93 t = obj.t; 68 t = obj.t;
94 end 69 end
95 70
96 function [vt,t] = getVt(obj) 71 function [vt,t] = getVt(obj)
97 % Calculate next time step to be able to do centered diff. 72 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
98 v_next = zeros(size(obj.v));
99 b = obj.G(obj.t) - obj.BB*obj.v - obj.CC*obj.v_prev;
100
101 y = obj.L\b(obj.p);
102 z = obj.U\y;
103 v_next(obj.q) = z;
104
105
106 vt = (v_next-obj.v_prev)/(2*obj.k);
107 t = obj.t; 73 t = obj.t;
108 end 74 end
109 75
110 % Calculate the conserved energy (Dm*v_n)^2_A + Im*v_n^2_B 76 % Calculate the conserved energy (Dm*v_n)^2_A + Im*v_n^2_B
111 function E = getEnergy(obj) 77 function E = getEnergy(obj)
112 v = obj.v; 78 v = obj.v;
113 vp = obj.v_prev; 79 vp = obj.v_prev;
114 vt = (obj.v - obj.v_prev)/obj.k; 80 vt = (obj.v - obj.v_prev)/obj.k;
115 81
116 E = vt'*obj.A*vt + 1/2*(v'*obj.B*v + vp'*obj.B*vp); 82 E = vt'*obj.A*vt + 1/2*(v'*obj.C*v + vp'*obj.C*vp);
117 end 83 end
118 84
119 function obj = step(obj) 85 function obj = step(obj)
120 b = obj.G(obj.t) - obj.BB*obj.v - obj.CC*obj.v_prev; 86 b = obj.G(obj.t) - obj.BB*obj.v - obj.CC*obj.v_prev;
121 obj.v_prev = obj.v; 87 obj.v_prev = obj.v;
132 obj.t = obj.t + obj.k; 98 obj.t = obj.t + obj.k;
133 obj.n = obj.n + 1; 99 obj.n = obj.n + 1;
134 end 100 end
135 end 101 end
136 end 102 end
137
138
139
140
141
142 %%% Derivation
143 % syms A B C G
144 % syms n k
145 % syms f1 f2
146
147 % v = symfun(sym('v(n)'),n);
148
149
150 % d = A/k^2 * (v(n+1) - 2*v(n) +v(n-1)) + B/2*(v(n+1)+v(n-1)) + C/(2*k)*(v(n+1) - v(n-1)) == G
151 % ic1 = v(0) == f1
152 % ic2 = A/k*(v(1)-f1) + k/2*(B*f1 + C*f2 - G) - f2 == 0
153
154 % c = collect(d, [v(n) v(n-1) v(n+1)]) % (-(2*A)/k^2)*v(n) + (B/2 + A/k^2 - C/(2*k))*v(n - 1) + (B/2 + A/k^2 + C/(2*k))*v(n + 1) == G
155 % syms AA BB CC
156 % % AA = B/2 + A/k^2 + C/(2*k)
157 % % BB = -(2*A)/k^2
158 % % CC = B/2 + A/k^2 - C/(2*k)
159 % s = subs(c, [B/2 + A/k^2 + C/(2*k), -(2*A)/k^2, B/2 + A/k^2 - C/(2*k)], [AA, BB, CC])
160
161
162 % ic2_a = collect(ic2, [v(1) f1 f2]) % (A/k)*v(1) + ((B*k)/2 - A/k)*f1 + ((C*k)/2 - 1)*f2 - (G*k)/2 == 0
163