Mercurial > repos > public > sbplib

comparison +time/Rungekutta4SecondOrder.m @ 886:8894e9c49e40 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default for latest changes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 15 Nov 2018 16:36:21 -0800
parents b5e5b195da1e ae905a11e32c
children 50d5a3843099
comparison
equal deleted inserted replaced
816:b5e5b195da1e 886:8894e9c49e40
13 n 13 n
14 end 14 end
15 15
16 16
17 methods 17 methods
18 % Solves u_tt = Du + Eu_t + S by
19 % Rewriting on first order form:
20 % w_t = M*w + C(t)
21 % where
22 % M = [
23 % 0, I;
24 % D, E;
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.
18 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) 33 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t)
19 obj.D = D; 34 obj.D = D;
20 obj.E = E; 35 obj.E = E;
21 obj.S = S; 36 obj.S = S;
22 obj.m = length(v0); 37 obj.m = length(v0);
23 obj.n = 0; 38 obj.n = 0;
24 39
25 I = speye(obj.m);
26 O = sparse(obj.m,obj.m);
27 obj.M = [O, I; D, E*I]; % Multiply with I to allow 0 as input.
28 40
29 if S == 0 41 if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle')
30 obj.C = zeros(2*obj.m,1); 42 default_arg('D', @(t)sparse(obj.m, obj.m));
43 default_arg('E', @(t)sparse(obj.m, obj.m));
44 default_arg('S', @(t)sparse(obj.m, 1) );
45
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 ];
31 else 65 else
32 obj.C = [zeros(obj.m,1), S]; 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);
33 end 88 end
34
35 obj.k = k;
36 obj.t = t0;
37 obj.w = [v0; v0t];
38
39 obj.F = @(w,t)(obj.M*w + obj.C);
40 end 89 end
41 90
42 function [v,t] = getV(obj) 91 function [v,t] = getV(obj)
43 v = obj.w(1:end/2); 92 v = obj.w(1:end/2);
44 t = obj.t; 93 t = obj.t;