Mercurial > repos > public > sbplib

comparison +time/Cdiff.m @ 977:c9009d5a3101

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor time.Cdiff and add getEnergyMethod
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 07 Jan 2019 15:56:35 +0100
parents 151ab2b5a686
children 1a30dbe99c7c
comparison
equal deleted inserted replaced
976:ea8fefc326b2 977:c9009d5a3101
1 classdef Cdiff < time.Timestepper 1 classdef Cdiff < time.Timestepper
2 properties 2 properties
3 D 3 A, B, C
4 E 4 AA, BB, CC
5 S 5 G
6 k 6 k
7 t 7 t
8 v 8 v
9 v_prev 9 v_prev
10 n 10 n
11 end 11 end
12 12
13 13
14 methods 14 methods
15 % Solves u_tt = Du + Eu_t + S 15 % Solves
16 % D, E, S can either all be constants or all be function handles, 16 % A*v_tt + B*v_t + C*v = G(t)
17 % They can also be omitted by setting them equal to the empty matrix. 17 % v(t0) = v0
18 % Cdiff(D, E, S, k, t0, n0, v, v_prev) 18 % v_t(t0) = v0t
19 function obj = Cdiff(D, E, S, k, t0, n0, v, v_prev) 19 % starting at time t0 with timestep k
20 m = length(v); 20 % Using
21 default_arg('E',sparse(m,m)); 21 % A*Dp*Dm*v_n + B*D0*v_n + C*v_n = G(t_n)
22 default_arg('S',sparse(m,1)); 22 function obj = Cdiff(A, B, C, G, v0, v0t, k, t0)
23 m = length(v0);
24 default_arg('A', speye(m));
25 default_arg('B', sparse(m,m));
26 default_arg('G', @(t) sparse(m,1));
27 default_arg('t0', 0);
23 28
24 obj.D = D; 29 obj.A = A;
25 obj.E = E; 30 obj.B = B;
26 obj.S = S; 31 obj.C = C;
32 obj.G = G;
27 33
34 % Rewrite as AA*v_(n+1) + BB*v_n + CC*v_(n-1) = G(t_n)
35 obj.AA = A/k^2 + B/2/k;
36 obj.BB = -2*A/k^2 + C;
37 obj.CC = A/k^2 - B/2/k;
28 38
29 obj.k = k; 39 obj.k = k;
30 obj.t = t0; 40 obj.v_prev = v0;
31 obj.n = n0; 41 obj.v = v0 + k*v0t;
32 obj.v = v; 42 obj.t = t0+k;
33 obj.v_prev = v_prev; 43 obj.n = 1;
34 end 44 end
35 45
36 function [v,t] = getV(obj) 46 function [v,t] = getV(obj)
37 v = obj.v; 47 v = obj.v;
38 t = obj.t; 48 t = obj.t;
41 function [vt,t] = getVt(obj) 51 function [vt,t] = getVt(obj)
42 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) 52 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
43 t = obj.t; 53 t = obj.t;
44 end 54 end
45 55
56 function E = getEnergy(obj)
57 v = obj.v;
58 vp = obj.v_prev;
59 vt = (obj.v - obj.v_prev)/obj.k;
60
61 E = vt'*obj.A*vt + v'*obj.C*vp;
62 end
63
46 function obj = step(obj) 64 function obj = step(obj)
47 [obj.v, obj.v_prev] = time.cdiff.cdiff(obj.v, obj.v_prev, obj.k, obj.D, obj.E, obj.S); 65 v_next = obj.AA\(-obj.BB*obj.v - obj.CC*obj.v_prev + obj.G(obj.t));
66
67 obj.v_prev = obj.v;
68 obj.v = v_next;
48 obj.t = obj.t + obj.k; 69 obj.t = obj.t + obj.k;
49 obj.n = obj.n + 1; 70 obj.n = obj.n + 1;
50 end 71 end
51 end 72 end
52 end 73 end