Mercurial > repos > public > sbplib

comparison +time/CdiffImplicit.m @ 425:e56dbd9e4196 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/beams
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 07 Feb 2017 16:09:02 +0100
parents 67ca8964f03c
children d5bce13ece23
comparison
equal deleted inserted replaced
423:a2cb0d4f4a02 425:e56dbd9e4196
1 classdef CdiffImplicit < time.Timestepper
2 properties
3 A, B, C, G
4 AA, BB, CC
5 k
6 t
7 v, v_prev
8 n
9
10 % LU factorization
11 L,U,p,q
12 end
13
14 methods
15 % Solves
16 % A*u_tt + B*u + C*v_t = G(t)
17 % u(t0) = f1
18 % u_t(t0) = f2
19 % starting at time t0 with timestep k
20 function obj = CdiffImplicit(A, B, C, G, f1, f2, k, t0)
21 default_arg('A', []);
22 default_arg('C', []);
23 default_arg('G', []);
24 default_arg('f1', 0);
25 default_arg('f2', 0);
26 default_arg('t0', 0);
27
28 m = size(B,1);
29
30 if isempty(A)
31 A = speye(m);
32 end
33
34 if isempty(C)
35 C = sparse(m,m);
36 end
37
38 if isempty(G)
39 G = @(t) sparse(m,1);
40 end
41
42 if isempty(f1)
43 f1 = sparse(m,1);
44 end
45
46 if isempty(f2)
47 f2 = sparse(m,1);
48 end
49
50 obj.A = A;
51 obj.B = B;
52 obj.C = C;
53 obj.G = G;
54
55 AA = 1/k^2*A + 1/2*B + 1/(2*k)*C;
56 BB = -2/k^2*A;
57 CC = 1/k^2*A + 1/2*B - 1/(2*k)*C;
58 % AA*v_next + BB*v + CC*v_prev == G(t_n)
59
60 obj.AA = AA;
61 obj.BB = BB;
62 obj.CC = CC;
63
64 v_prev = f1;
65 I = speye(m);
66 v = (1/k^2*A)\((1/k^2*A - 1/2*B)*f1 + (1/k*I - 1/2*C)*f2 + 1/2*G(0));
67
68 if ~issparse(A) || ~issparse(B) || ~issparse(C)
69 error('LU factorization with full pivoting only works for sparse matrices.')
70 end
71
72 [L,U,p,q] = lu(AA,'vector');
73
74 obj.L = L;
75 obj.U = U;
76 obj.p = p;
77 obj.q = q;
78
79
80 obj.k = k;
81 obj.t = t0+k;
82 obj.n = 1;
83 obj.v = v;
84 obj.v_prev = v_prev;
85 end
86
87 function [v,t] = getV(obj)
88 v = obj.v;
89 t = obj.t;
90 end
91
92 function [vt,t] = getVt(obj)
93 % Calculate next time step to be able to do centered diff.
94 v_next = zeros(size(obj.v));
95 b = obj.G(obj.t) - obj.BB*obj.v - obj.CC*obj.v_prev;
96
97 y = obj.L\b(obj.p);
98 z = obj.U\y;
99 v_next(obj.q) = z;
100
101
102 vt = (v_next-obj.v_prev)/(2*obj.k);
103 t = obj.t;
104 end
105
106 function obj = step(obj)
107 b = obj.G(obj.t) - obj.BB*obj.v - obj.CC*obj.v_prev;
108 obj.v_prev = obj.v;
109
110 % % Backslash
111 % obj.v = obj.AA\b;
112
113 % LU with column pivot
114 y = obj.L\b(obj.p);
115 z = obj.U\y;
116 obj.v(obj.q) = z;
117
118 % Update time
119 obj.t = obj.t + obj.k;
120 obj.n = obj.n + 1;
121 end
122 end
123 end
124
125
126
127
128
129 %%% Derivation
130 % syms A B C G
131 % syms n k
132 % syms f1 f2
133
134 % v = symfun(sym('v(n)'),n);
135
136
137 % 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
138 % ic1 = v(0) == f1
139 % ic2 = A/k*(v(1)-f1) + k/2*(B*f1 + C*f2 - G) - f2 == 0
140
141 % 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
142 % syms AA BB CC
143 % % AA = B/2 + A/k^2 + C/(2*k)
144 % % BB = -(2*A)/k^2
145 % % CC = B/2 + A/k^2 - C/(2*k)
146 % 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])
147
148
149 % 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
150