Mercurial > repos > public > sbplib

changeset 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 4d472d020ccf
children 3a286d2d1939
files +time/Rungekutta4SecondOrder.m
diffstat 1 files changed, 30 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
diff -r 4d472d020ccf -r 74eec7e69b63 +time/Rungekutta4SecondOrder.m
--- a/+time/Rungekutta4SecondOrder.m	Wed Feb 12 13:27:27 2020 +0100
+++ b/+time/Rungekutta4SecondOrder.m	Fri Mar 13 12:07:00 2020 +0100
@@ -36,9 +36,32 @@
             obj.S = S;
             obj.m = length(v0);
             obj.n = 0;
-
+            
+            % Construct RHS system on first order form. 
+            % 3 different setups are handeled:
+            % If none of D, E, S are time-dependent (e.g function_handles)
+            % then a complete matrix of the RHS is constructed.
+            % If the source term S is time-depdentent, then the first order system matrix M
+            % is formed via D and E, while S is kept as a function handle
+            % If D or E are time-dependent, then all terns are treated as function handles.
+            if ~isa(D, 'function_handle') && ~isa(E, 'function_handle') && isa(S, 'function_handle')
+                default_arg('D', sparse(obj.m, obj.m));
+                default_arg('E', sparse(obj.m, obj.m));
+                default_arg('S', @(t)sparse(obj.m, 1));
+                
+                I = speye(obj.m);
+                O = sparse(obj.m,obj.m);
 
-            if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle')
+                obj.M = [
+                    O, I;
+                    D, E;
+                ];
+
+                obj.k = k;
+                obj.t = t0;
+                obj.w = [v0; v0t];
+                obj.F = @(w,t)obj.rhsTimeDepS(w,t,obj.M,S,obj.m);
+            elseif isa(D, 'function_handle') || isa(E, 'function_handle')
                 default_arg('D', @(t)sparse(obj.m, obj.m));
                 default_arg('E', @(t)sparse(obj.m, obj.m));
                 default_arg('S', @(t)sparse(obj.m, 1)    );
@@ -63,7 +86,6 @@
                     D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t);
                 ];
             else
-
                 default_arg('D', sparse(obj.m, obj.m));
                 default_arg('E', sparse(obj.m, obj.m));
                 default_arg('S', sparse(obj.m, 1)    );
@@ -110,6 +132,11 @@
         function k = getTimeStep(lambda)
             k = rk4.get_rk4_time_step(lambda);
         end
+        
+        function F = rhsTimeDepS(w,t,M,S,m)
+            F = M*w;
+            F(m+1:end) = F(m+1:end) + S(t);
+        end
     end
 
 end
\ No newline at end of file