Mercurial > repos > public > sbplib

comparison +time/Rungekutta4SecondOrder.m @ 425:e56dbd9e4196 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/beams
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 07 Feb 2017 16:09:02 +0100
parents 7b5ef8b89268
children ae905a11e32c
comparison
equal deleted inserted replaced
423:a2cb0d4f4a02 425:e56dbd9e4196
13 n 13 n
14 end 14 end
15 15
16 16
17 methods 17 methods
18 % Solves u_tt = Du + Eu_t + S
19 % D, E, S can either all be constants or all be function handles,
20 % They can also be omitted by setting them equal to the empty matrix.
18 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) 21 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t)
19 obj.D = D; 22 obj.D = D;
20 obj.E = E; 23 obj.E = E;
21 obj.S = S; 24 obj.S = S;
22 obj.m = length(v0); 25 obj.m = length(v0);
23 obj.n = 0; 26 obj.n = 0;
24 27
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 28
29 if S == 0 29 if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle')
30 obj.C = zeros(2*obj.m,1); 30 default_arg('D', @(t)sparse(obj.m, obj.m));
31 default_arg('E', @(t)sparse(obj.m, obj.m));
32 default_arg('S', @(t)sparse(obj.m, 1) );
33
34 if ~isa(D, 'function_handle')
35 D = @(t)D;
36 end
37 if ~isa(E, 'function_handle')
38 E = @(t)E;
39 end
40 if ~isa(S, 'function_handle')
41 S = @(t)S;
42 end
43
44 I = speye(obj.m);
45 O = sparse(obj.m,obj.m);
46
47 obj.M = @(t)[
48 O, I;
49 D(t), E(t);
50 ];
51 obj.C = @(t)[
52 zeros(obj.m,1);
53 S(t);
54 ];
55
56 obj.k = k;
57 obj.t = t0;
58 obj.w = [v0; v0t];
59
60 obj.F = @(w,t)(obj.M(t)*w + obj.C(t));
31 else 61 else
32 obj.C = [zeros(obj.m,1), S]; 62
63 default_arg('D', sparse(obj.m, obj.m));
64 default_arg('E', sparse(obj.m, obj.m));
65 default_arg('S', sparse(obj.m, 1) );
66
67 I = speye(obj.m);
68 O = sparse(obj.m,obj.m);
69
70 obj.M = [
71 O, I;
72 D, E;
73 ];
74 obj.C = [
75 zeros(obj.m,1);
76 S;
77 ];
78
79 obj.k = k;
80 obj.t = t0;
81 obj.w = [v0; v0t];
82
83 obj.F = @(w,t)(obj.M*w + obj.C);
33 end 84 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 85 end
41 86
42 function [v,t] = getV(obj) 87 function [v,t] = getV(obj)
43 v = obj.w(1:end/2); 88 v = obj.w(1:end/2);
44 t = obj.t; 89 t = obj.t;