Mercurial > repos > public > sbplib

changeset 984:0585a2ee7ee7 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Inline the rk.rungekutta_4 function.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 08 Jan 2019 12:19:33 +0100
parents b89379fb0814
children ad6de007e7f6
files +time/+rk/rungekutta_4.m +time/Rk4SecondOrderNonlin.m +time/Rungekutta.m +time/Rungekutta4.m +time/Rungekutta4SecondOrder.m
diffstat 5 files changed, 28 insertions(+), 20 deletions(-) [+]
line wrap: on
line diff
diff -r b89379fb0814 -r 0585a2ee7ee7 +time/+rk/rungekutta_4.m
--- a/+time/+rk/rungekutta_4.m	Tue Jan 08 12:07:57 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,10 +0,0 @@
-% Takes one time step of size k using the rungekutta method
-% starting from v_0 and where the function F(v,t) gives the
-% time derivatives.
-function v = rungekutta_4(v, t , k, F)
-    k1 = F(v         ,t      );
-    k2 = F(v+0.5*k*k1,t+0.5*k);
-    k3 = F(v+0.5*k*k2,t+0.5*k);
-    k4 = F(v+    k*k3,t+    k);
-    v = v + (1/6)*(k1+2*(k2+k3)+k4)*k;
-end
\ No newline at end of file
diff -r b89379fb0814 -r 0585a2ee7ee7 +time/Rk4SecondOrderNonlin.m
--- a/+time/Rk4SecondOrderNonlin.m	Tue Jan 08 12:07:57 2019 +0100
+++ b/+time/Rk4SecondOrderNonlin.m	Tue Jan 08 12:19:33 2019 +0100
@@ -61,7 +61,15 @@
         end
 
         function obj = step(obj)
-            obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
+            w = obj.w;
+            k = obj.k;
+
+            k1 = obj.F(t, w);
+            k2 = obj.F(t + 0.5*k, w + 0.5*k*k1);
+            k3 = obj.F(t + 0.5*k, w + 0.5*k*k2);
+            k4 = obj.F(t + k, w + k*k3);
+
+            obj.w = w + k*(1/6)*(k1+2*(k2+k3)+k4);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end
diff -r b89379fb0814 -r 0585a2ee7ee7 +time/Rungekutta.m
--- a/+time/Rungekutta.m	Tue Jan 08 12:07:57 2019 +0100
+++ b/+time/Rungekutta.m	Tue Jan 08 12:19:33 2019 +0100
@@ -30,13 +30,7 @@
             obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
 
             if isempty(discreteData)
-                % TODO: method "rk4" is also implemented in the butcher tableau, but the rungekutta_4.m implementation
-                % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it.
-                if (method == "rk4")
-                    obj.scheme = @(v,t,n) time.rk.rungekutta_4(v ,t, dt, F);
-                else
-                    obj.scheme = @(v,t,n) time.rk.rungekutta(v, t, dt, F, obj.coeffs);
-                end
+                obj.scheme = @(v,t,n) time.rk.rungekutta(v, t, dt, F, obj.coeffs);
             else
                 obj.scheme = @(v,t,n) time.rk.rungekuttaDiscreteData(v, t, dt, F, obj.coeffs, discreteData, n);
             end
diff -r b89379fb0814 -r 0585a2ee7ee7 +time/Rungekutta4.m
--- a/+time/Rungekutta4.m	Tue Jan 08 12:07:57 2019 +0100
+++ b/+time/Rungekutta4.m	Tue Jan 08 12:19:33 2019 +0100
@@ -26,7 +26,15 @@
         end
 
         function obj = step(obj)
-            obj.v = time.rk.rungekutta_4(obj.v, obj.t, obj.dt, obj.F);
+            v = obj.v;
+            dt = obj.dt;
+
+            k1 = obj.F(t, v);
+            k2 = obj.F(t + 0.5*dt, v + 0.5*dt*k1);
+            k3 = obj.F(t + 0.5*dt, v + 0.5*dt*k2);
+            k4 = obj.F(t + dt, v + dt*k3);
+
+            obj.v = v + dt*(1/6)*(k1+2*(k2+k3)+k4);
             obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
diff -r b89379fb0814 -r 0585a2ee7ee7 +time/Rungekutta4SecondOrder.m
--- a/+time/Rungekutta4SecondOrder.m	Tue Jan 08 12:07:57 2019 +0100
+++ b/+time/Rungekutta4SecondOrder.m	Tue Jan 08 12:19:33 2019 +0100
@@ -99,7 +99,15 @@
         end
 
         function obj = step(obj)
-            obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
+            w = obj.w;
+            k = obj.k;
+
+            k1 = obj.F(t, w);
+            k2 = obj.F(t + 0.5*k, w + 0.5*k*k1);
+            k3 = obj.F(t + 0.5*k, w + 0.5*k*k2);
+            k4 = obj.F(t + k, w + k*k3);
+
+            obj.w = w + k*(1/6)*(k1+2*(k2+k3)+k4);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end