Mercurial > repos > public > sbplib

comparison +time/Cdiff.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 7a5e770974ed
children
comparison
equal deleted inserted replaced
1112:835c8fa456ec 1113:47e86b5270ad
1 classdef Cdiff < time.Timestepper 1 classdef Cdiff < 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 end 10 end
11 11
13 methods 13 methods
14 % Solves 14 % Solves
15 % A*v_tt + B*v_t + C*v = G(t) 15 % A*v_tt + B*v_t + C*v = G(t)
16 % v(t0) = v0 16 % v(t0) = v0
17 % v_t(t0) = v0t 17 % v_t(t0) = v0t
18 % starting at time t0 with timestep k 18 % starting at time t0 with timestep dt
19 % Using 19 % Using
20 % A*Dp*Dm*v_n + B*D0*v_n + C*v_n = G(t_n) 20 % A*Dp*Dm*v_n + B*D0*v_n + C*v_n = G(t_n)
21 function obj = Cdiff(A, B, C, G, v0, v0t, k, t0) 21 function obj = Cdiff(A, B, C, G, v0, v0t, dt, t0)
22 m = length(v0); 22 m = length(v0);
23 default_arg('A', speye(m)); 23 default_arg('A', speye(m));
24 default_arg('B', sparse(m,m)); 24 default_arg('B', sparse(m,m));
25 default_arg('G', @(t) sparse(m,1)); 25 default_arg('G', @(t) sparse(m,1));
26 default_arg('t0', 0); 26 default_arg('t0', 0);
29 obj.B = B; 29 obj.B = B;
30 obj.C = C; 30 obj.C = C;
31 obj.G = G; 31 obj.G = G;
32 32
33 % Rewrite as AA*v_(n+1) + BB*v_n + CC*v_(n-1) = G(t_n) 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); 34 obj.AA = A/dt^2 + B/(2*dt);
35 obj.BB = -2*A/k^2 + C; 35 obj.BB = -2*A/dt^2 + C;
36 obj.CC = A/k^2 - B/(2*k); 36 obj.CC = A/dt^2 - B/(2*dt);
37 37
38 obj.k = k; 38 obj.dt = dt;
39 obj.v_prev = v0; 39 obj.v_prev = v0;
40 obj.v = v0 + k*v0t; 40 obj.v = v0 + dt*v0t;
41 obj.t = t0+k; 41 obj.t = t0+dt;
42 obj.n = 1; 42 obj.n = 1;
43 end 43 end
44 44
45 function [v,t] = getV(obj) 45 function [v,t] = getV(obj)
46 v = obj.v; 46 v = obj.v;
47 t = obj.t; 47 t = obj.t;
48 end 48 end
49 49
50 function [vt,t] = getVt(obj) 50 function [vt,t] = getVt(obj)
51 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) 51 vt = (obj.v-obj.v_prev)/obj.dt; % Could be improved using u_tt = f(u))
52 t = obj.t; 52 t = obj.t;
53 end 53 end
54 54
55 function E = getEnergy(obj) 55 function E = getEnergy(obj)
56 v = obj.v; 56 v = obj.v;
57 vp = obj.v_prev; 57 vp = obj.v_prev;
58 vt = (obj.v - obj.v_prev)/obj.k; 58 vt = (obj.v - obj.v_prev)/obj.dt;
59 59
60 E = vt'*obj.A*vt + v'*obj.C*vp; 60 E = vt'*obj.A*vt + v'*obj.C*vp;
61 end 61 end
62 62
63 function obj = step(obj) 63 function obj = step(obj)
64 v_next = obj.AA\(-obj.BB*obj.v - obj.CC*obj.v_prev + obj.G(obj.t)); 64 v_next = obj.AA\(-obj.BB*obj.v - obj.CC*obj.v_prev + obj.G(obj.t));
65 65
66 obj.v_prev = obj.v; 66 obj.v_prev = obj.v;
67 obj.v = v_next; 67 obj.v = v_next;
68 obj.t = obj.t + obj.k; 68 obj.t = obj.t + obj.dt;
69 obj.n = obj.n + 1; 69 obj.n = obj.n + 1;
70 end 70 end
71 end 71 end
72 end 72 end