Mercurial > repos > public > sbplib

comparison +time/Rungekutta4SecondOrder.m @ 1258:74eec7e69b63 feature/FMMlabb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 13 Mar 2020 12:07:00 +0100
parents ae905a11e32c
children
comparison
equal deleted inserted replaced
1257:4d472d020ccf 1258:74eec7e69b63
34 obj.D = D; 34 obj.D = D;
35 obj.E = E; 35 obj.E = E;
36 obj.S = S; 36 obj.S = S;
37 obj.m = length(v0); 37 obj.m = length(v0);
38 obj.n = 0; 38 obj.n = 0;
39
40 % Construct RHS system on first order form.
41 % 3 different setups are handeled:
42 % If none of D, E, S are time-dependent (e.g function_handles)
43 % then a complete matrix of the RHS is constructed.
44 % If the source term S is time-depdentent, then the first order system matrix M
45 % is formed via D and E, while S is kept as a function handle
46 % If D or E are time-dependent, then all terns are treated as function handles.
47 if ~isa(D, 'function_handle') && ~isa(E, 'function_handle') && isa(S, 'function_handle')
48 default_arg('D', sparse(obj.m, obj.m));
49 default_arg('E', sparse(obj.m, obj.m));
50 default_arg('S', @(t)sparse(obj.m, 1));
51
52 I = speye(obj.m);
53 O = sparse(obj.m,obj.m);
39 54
55 obj.M = [
56 O, I;
57 D, E;
58 ];
40 59
41 if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle') 60 obj.k = k;
61 obj.t = t0;
62 obj.w = [v0; v0t];
63 obj.F = @(w,t)obj.rhsTimeDepS(w,t,obj.M,S,obj.m);
64 elseif isa(D, 'function_handle') || isa(E, 'function_handle')
42 default_arg('D', @(t)sparse(obj.m, obj.m)); 65 default_arg('D', @(t)sparse(obj.m, obj.m));
43 default_arg('E', @(t)sparse(obj.m, obj.m)); 66 default_arg('E', @(t)sparse(obj.m, obj.m));
44 default_arg('S', @(t)sparse(obj.m, 1) ); 67 default_arg('S', @(t)sparse(obj.m, 1) );
45 68
46 if ~isa(D, 'function_handle') 69 if ~isa(D, 'function_handle')
61 obj.F = @(w,t)[ 84 obj.F = @(w,t)[
62 w(obj.m+1:end); 85 w(obj.m+1:end);
63 D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t); 86 D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t);
64 ]; 87 ];
65 else 88 else
66
67 default_arg('D', sparse(obj.m, obj.m)); 89 default_arg('D', sparse(obj.m, obj.m));
68 default_arg('E', sparse(obj.m, obj.m)); 90 default_arg('E', sparse(obj.m, obj.m));
69 default_arg('S', sparse(obj.m, 1) ); 91 default_arg('S', sparse(obj.m, 1) );
70 92
71 I = speye(obj.m); 93 I = speye(obj.m);
108 130
109 methods (Static) 131 methods (Static)
110 function k = getTimeStep(lambda) 132 function k = getTimeStep(lambda)
111 k = rk4.get_rk4_time_step(lambda); 133 k = rk4.get_rk4_time_step(lambda);
112 end 134 end
135
136 function F = rhsTimeDepS(w,t,M,S,m)
137 F = M*w;
138 F(m+1:end) = F(m+1:end) + S(t);
139 end
113 end 140 end
114 141
115 end 142 end