Mercurial > repos > public > sbplib

changeset 991:a99f00896b8e feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Rungekutta4SecondOrder native second order
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 09 Jan 2019 10:17:00 +0100
parents 1066bb31bc95
children bbd165cc585c
files +time/Rungekutta4.m +time/Rungekutta4SecondOrder.m +time/Timestepper.m
diffstat 3 files changed, 18 insertions(+), 24 deletions(-) [+]
line wrap: on
line diff
--- a/+time/Rungekutta4.m	Wed Jan 09 09:09:15 2019 +0100
+++ b/+time/Rungekutta4.m	Wed Jan 09 10:17:00 2019 +0100
@@ -26,6 +26,7 @@
         end
 
         function obj = step(obj)
+            t = obj.t;
             v = obj.v;
             dt = obj.dt;
 
--- a/+time/Rungekutta4SecondOrder.m	Wed Jan 09 09:09:15 2019 +0100
+++ b/+time/Rungekutta4SecondOrder.m	Wed Jan 09 10:17:00 2019 +0100
@@ -2,11 +2,8 @@
     properties
         F
         dt
-        t
-        n
-
-        G % RHS for rewritten system
-        w
+        t, n
+        v, v_t
     end
 
 
@@ -20,45 +17,41 @@
             obj.t = t0;
             obj.n = 0;
 
-            m = length(v0);
-            obj.w = [v0; v0t];
-            obj.G = @(t, w)[
-                w(m+1:end);
-                F(t,w(1:m), w(m+1:end));
-            ];
+            obj.v = v0;
+            obj.v_t = v0t;
         end
 
         function [v,t] = getV(obj)
-            v = obj.w(1:end/2);
+            v = obj.v
             t = obj.t;
         end
 
         function [vt,t] = getVt(obj)
-            vt = obj.w(end/2+1:end);
+            vt = obj.v_t;
             t = obj.t;
         end
 
         function obj = step(obj)
-            % TODO: Avoid extra storage
-            w = obj.w;
+            t = obj.t;
+            v = obj.v;
+            v_t = obj.v_t;
             dt = obj.dt;
 
-            k1 = obj.G(t, w);
-            k2 = obj.G(t + 0.5*dt, w + 0.5*dt*k1);
-            k3 = obj.G(t + 0.5*dt, w + 0.5*dt*k2);
-            k4 = obj.G(t + dt, w + dt*k3);
+            k1 = obj.F(t,          v,                            v_t);
+            k2 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t,               v_t + 1/2*dt*k1);
+            k3 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t + 1/4*dt^2*k1, v_t + 1/2*dt*k2);
+            k4 = obj.F(t + dt,     v + dt*v_t + 1/2*dt^2*k2,     v_t + dt*k3);
 
-            obj.w = w + dt*(1/6)*(k1+2*(k2+k3)+k4);
+            obj.v = v + dt*v_t + dt^2*(1/6)*(k1 + k2 + k3);
+            obj.v_t = v_t + dt*(1/6)*(k1 + 2*k2 + 2*k3 + k4);
             obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
 
-
     methods (Static)
         function k = getTimeStep(lambda)
             k = rk4.get_rk4_time_step(lambda);
         end
     end
-
-end
\ No newline at end of file
+end
--- a/+time/Timestepper.m	Wed Jan 09 09:09:15 2019 +0100
+++ b/+time/Timestepper.m	Wed Jan 09 10:17:00 2019 +0100
@@ -1,7 +1,7 @@
 classdef Timestepper < handle
     properties (Abstract)
         t
-        k
+        k % TBD: should this be a method instead?
         n
     end