Mercurial > repos > public > sbplib

comparison +time/CdiffImplicit.m @ 1113:47e86b5270ad feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change name of property k to dt in time.Timestepper
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 10 Apr 2019 22:40:55 +0200
parents 1a30dbe99c7c
children
comparison
equal deleted inserted replaced
1112:835c8fa456ec 1113:47e86b5270ad
1 classdef CdiffImplicit < time.Timestepper 1 classdef CdiffImplicit < time.Timestepper
2 properties 2 properties
3 A, B, C 3 A, B, C
4 AA, BB, CC 4 AA, BB, CC
5 G 5 G
6 k 6 dt
7 t 7 t
8 v, v_prev 8 v, v_prev
9 n 9 n
10 10
11 % LU factorization 11 % LU factorization
18 % v(t0) = v0 18 % v(t0) = v0
19 % v_t(t0) = v0t 19 % v_t(t0) = v0t
20 % starting at time t0 with timestep 20 % starting at time t0 with timestep
21 % Using 21 % Using
22 % A*Dp*Dm*v_n + B*D0*v_n + C*I0*v_n = G(t_n) 22 % A*Dp*Dm*v_n + B*D0*v_n + C*I0*v_n = G(t_n)
23 function obj = CdiffImplicit(A, B, C, G, v0, v0t, k, t0) 23 function obj = CdiffImplicit(A, B, C, G, v0, v0t, dt, t0)
24 m = length(v0); 24 m = length(v0);
25 default_arg('A', speye(m)); 25 default_arg('A', speye(m));
26 default_arg('B', sparse(m,m)); 26 default_arg('B', sparse(m,m));
27 default_arg('G', @(t) sparse(m,1)); 27 default_arg('G', @(t) sparse(m,1));
28 default_arg('t0', 0); 28 default_arg('t0', 0);
31 obj.B = B; 31 obj.B = B;
32 obj.C = C; 32 obj.C = C;
33 obj.G = G; 33 obj.G = G;
34 34
35 % Rewrite as AA*v_(n+1) + BB*v_n + CC*v_(n-1) = G(t_n) 35 % Rewrite as AA*v_(n+1) + BB*v_n + CC*v_(n-1) = G(t_n)
36 AA = A/k^2 + B/(2*k) + C/2; 36 AA = A/dt^2 + B/(2*dt) + C/2;
37 BB = -2*A/k^2; 37 BB = -2*A/dt^2;
38 CC = A/k^2 - B/(2*k) + C/2; 38 CC = A/dt^2 - B/(2*dt) + C/2;
39 39
40 obj.AA = AA; 40 obj.AA = AA;
41 obj.BB = BB; 41 obj.BB = BB;
42 obj.CC = CC; 42 obj.CC = CC;
43 43
44 v_prev = v0; 44 v_prev = v0;
45 I = speye(m); 45 I = speye(m);
46 v = v0 + k*v0t; 46 v = v0 + dt*v0t;
47 47
48 if ~issparse(A) || ~issparse(B) || ~issparse(C) 48 if ~issparse(A) || ~issparse(B) || ~issparse(C)
49 error('LU factorization with full pivoting only works for sparse matrices.') 49 error('LU factorization with full pivoting only works for sparse matrices.')
50 end 50 end
51 51
54 obj.L = L; 54 obj.L = L;
55 obj.U = U; 55 obj.U = U;
56 obj.p = p; 56 obj.p = p;
57 obj.q = q; 57 obj.q = q;
58 58
59 obj.k = k; 59 obj.dt = dt;
60 obj.t = t0+k; 60 obj.t = t0+dt;
61 obj.n = 1; 61 obj.n = 1;
62 obj.v = v; 62 obj.v = v;
63 obj.v_prev = v_prev; 63 obj.v_prev = v_prev;
64 end 64 end
65 65
67 v = obj.v; 67 v = obj.v;
68 t = obj.t; 68 t = obj.t;
69 end 69 end
70 70
71 function [vt,t] = getVt(obj) 71 function [vt,t] = getVt(obj)
72 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) 72 vt = (obj.v-obj.v_prev)/obj.dt; % Could be improved using u_tt = f(u))
73 t = obj.t; 73 t = obj.t;
74 end 74 end
75 75
76 % 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
77 function E = getEnergy(obj) 77 function E = getEnergy(obj)
78 v = obj.v; 78 v = obj.v;
79 vp = obj.v_prev; 79 vp = obj.v_prev;
80 vt = (obj.v - obj.v_prev)/obj.k; 80 vt = (obj.v - obj.v_prev)/obj.dt;
81 81
82 E = vt'*obj.A*vt + 1/2*(v'*obj.C*v + vp'*obj.C*vp); 82 E = vt'*obj.A*vt + 1/2*(v'*obj.C*v + vp'*obj.C*vp);
83 end 83 end
84 84
85 function obj = step(obj) 85 function obj = step(obj)
93 y = obj.L\b(obj.p); 93 y = obj.L\b(obj.p);
94 z = obj.U\y; 94 z = obj.U\y;
95 obj.v(obj.q) = z; 95 obj.v(obj.q) = z;
96 96
97 % Update time 97 % Update time
98 obj.t = obj.t + obj.k; 98 obj.t = obj.t + obj.dt;
99 obj.n = obj.n + 1; 99 obj.n = obj.n + 1;
100 end 100 end
101 end 101 end
102 end 102 end