Mercurial > repos > public > sbplib

comparison +time/Rungekutta4SecondOrder.m @ 707:0de70ec8bf60 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge with feature/optim
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 10 Nov 2017 14:22:56 +0100
parents ae905a11e32c
children c6fcee3fcf1b 8894e9c49e40 74eec7e69b63
comparison
equal deleted inserted replaced
696:7c16b5af8d98 707:0de70ec8bf60
13 n 13 n
14 end 14 end
15 15
16 16
17 methods 17 methods
18 % Solves u_tt = Du + Eu_t + S 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 % ]
19 % D, E, S can either all be constants or all be function handles, 31 % 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. 32 % They can also be omitted by setting them equal to the empty matrix.
21 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) 33 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t)
22 obj.D = D; 34 obj.D = D;
23 obj.E = E; 35 obj.E = E;
39 end 51 end
40 if ~isa(S, 'function_handle') 52 if ~isa(S, 'function_handle')
41 S = @(t)S; 53 S = @(t)S;
42 end 54 end
43 55
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; 56 obj.k = k;
57 obj.t = t0; 57 obj.t = t0;
58 obj.w = [v0; v0t]; 58 obj.w = [v0; v0t];
59 59
60 obj.F = @(w,t)(obj.M(t)*w + obj.C(t)); 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 ];
61 else 65 else
62 66
63 default_arg('D', sparse(obj.m, obj.m)); 67 default_arg('D', sparse(obj.m, obj.m));
64 default_arg('E', sparse(obj.m, obj.m)); 68 default_arg('E', sparse(obj.m, obj.m));
65 default_arg('S', sparse(obj.m, 1) ); 69 default_arg('S', sparse(obj.m, 1) );