Mercurial > repos > public > sbplib

comparison +time/Rungekutta4SecondOrder.m @ 985:ad6de007e7f6 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Convert Rungekutta4SecondOrder to take F(t,v,v_t) instead of matrices
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 08 Jan 2019 12:34:29 +0100
parents 0585a2ee7ee7
children a99f00896b8e
comparison
equal deleted inserted replaced
984:0585a2ee7ee7 985:ad6de007e7f6
1 classdef Rungekutta4SecondOrder < time.Timestepper 1 classdef Rungekutta4SecondOrder < time.Timestepper
2 properties 2 properties
3 F 3 F
4 k 4 dt
5 t 5 t
6 n
7
8 G % RHS for rewritten system
6 w 9 w
7 m
8 D
9 E
10 S
11 M
12 C
13 n
14 end 10 end
15 11
16 12
17 methods 13 methods
18 % Solves u_tt = Du + Eu_t + S by 14 % Create a time stepper for
19 % Rewriting on first order form: 15 % v_tt = F(t,v,v_t), v(t0) = v0, v_t(t0) = v0t
20 % w_t = M*w + C(t) 16 % with step size dt, by rewriting on first order form
21 % where 17 function obj = Rungekutta4SecondOrder(F, dt, t0, v0, v0t)
22 % M = [ 18 obj.F = F;
23 % 0, I; 19 obj.dt = dt;
24 % D, E; 20 obj.t = t0;
25 % ]
26 % and
27 % C(t) = [
28 % 0;
29 % S(t)
30 % ]
31 % D, E, S can either all be constants or all be function handles,
32 % They can also be omitted by setting them equal to the empty matrix.
33 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t)
34 obj.D = D;
35 obj.E = E;
36 obj.S = S;
37 obj.m = length(v0);
38 obj.n = 0; 21 obj.n = 0;
39 22
40 23 m = length(v0);
41 if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle') 24 obj.w = [v0; v0t];
42 default_arg('D', @(t)sparse(obj.m, obj.m)); 25 obj.G = @(t, w)[
43 default_arg('E', @(t)sparse(obj.m, obj.m)); 26 w(m+1:end);
44 default_arg('S', @(t)sparse(obj.m, 1) ); 27 F(t,w(1:m), w(m+1:end));
45 28 ];
46 if ~isa(D, 'function_handle')
47 D = @(t)D;
48 end
49 if ~isa(E, 'function_handle')
50 E = @(t)E;
51 end
52 if ~isa(S, 'function_handle')
53 S = @(t)S;
54 end
55
56 obj.k = k;
57 obj.t = t0;
58 obj.w = [v0; v0t];
59
60 % Avoid matrix formulation because it is VERY slow
61 obj.F = @(w,t)[
62 w(obj.m+1:end);
63 D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t);
64 ];
65 else
66
67 default_arg('D', sparse(obj.m, obj.m));
68 default_arg('E', sparse(obj.m, obj.m));
69 default_arg('S', sparse(obj.m, 1) );
70
71 I = speye(obj.m);
72 O = sparse(obj.m,obj.m);
73
74 obj.M = [
75 O, I;
76 D, E;
77 ];
78 obj.C = [
79 zeros(obj.m,1);
80 S;
81 ];
82
83 obj.k = k;
84 obj.t = t0;
85 obj.w = [v0; v0t];
86
87 obj.F = @(w,t)(obj.M*w + obj.C);
88 end
89 end 29 end
90 30
91 function [v,t] = getV(obj) 31 function [v,t] = getV(obj)
92 v = obj.w(1:end/2); 32 v = obj.w(1:end/2);
93 t = obj.t; 33 t = obj.t;
97 vt = obj.w(end/2+1:end); 37 vt = obj.w(end/2+1:end);
98 t = obj.t; 38 t = obj.t;
99 end 39 end
100 40
101 function obj = step(obj) 41 function obj = step(obj)
42 % TODO: Avoid extra storage
102 w = obj.w; 43 w = obj.w;
103 k = obj.k; 44 dt = obj.dt;
104 45
105 k1 = obj.F(t, w); 46 k1 = obj.G(t, w);
106 k2 = obj.F(t + 0.5*k, w + 0.5*k*k1); 47 k2 = obj.G(t + 0.5*dt, w + 0.5*dt*k1);
107 k3 = obj.F(t + 0.5*k, w + 0.5*k*k2); 48 k3 = obj.G(t + 0.5*dt, w + 0.5*dt*k2);
108 k4 = obj.F(t + k, w + k*k3); 49 k4 = obj.G(t + dt, w + dt*k3);
109 50
110 obj.w = w + k*(1/6)*(k1+2*(k2+k3)+k4); 51 obj.w = w + dt*(1/6)*(k1+2*(k2+k3)+k4);
111 obj.t = obj.t + obj.k; 52 obj.t = obj.t + obj.dt;
112 obj.n = obj.n + 1; 53 obj.n = obj.n + 1;
113 end 54 end
114 end 55 end
115 56
116 57