Mercurial > repos > public > sbplib

comparison +time/Cdiff.m @ 979:7a5e770974ed feature/timesteppers

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