diff +time/CdiffImplicit.m @ 1113:47e86b5270ad feature/timesteppers

Change name of property k to dt in time.Timestepper
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 10 Apr 2019 22:40:55 +0200
parents 1a30dbe99c7c
children
line wrap: on
line diff
--- a/+time/CdiffImplicit.m	Wed Apr 10 22:22:46 2019 +0200
+++ b/+time/CdiffImplicit.m	Wed Apr 10 22:40:55 2019 +0200
@@ -3,7 +3,7 @@
         A, B, C
         AA, BB, CC
         G
-        k
+        dt
         t
         v, v_prev
         n
@@ -20,7 +20,7 @@
         % starting at time t0 with timestep
         % Using
         % A*Dp*Dm*v_n + B*D0*v_n + C*I0*v_n = G(t_n)
-        function obj = CdiffImplicit(A, B, C, G, v0, v0t, k, t0)
+        function obj = CdiffImplicit(A, B, C, G, v0, v0t, dt, t0)
             m = length(v0);
             default_arg('A', speye(m));
             default_arg('B', sparse(m,m));
@@ -33,9 +33,9 @@
             obj.G = G;
 
             % Rewrite as AA*v_(n+1) + BB*v_n + CC*v_(n-1) = G(t_n)
-            AA =    A/k^2 + B/(2*k) + C/2;
-            BB = -2*A/k^2;
-            CC =    A/k^2 - B/(2*k) + C/2;
+            AA =    A/dt^2 + B/(2*dt) + C/2;
+            BB = -2*A/dt^2;
+            CC =    A/dt^2 - B/(2*dt) + C/2;
 
             obj.AA = AA;
             obj.BB = BB;
@@ -43,7 +43,7 @@
 
             v_prev = v0;
             I = speye(m);
-            v = v0 + k*v0t;
+            v = v0 + dt*v0t;
 
             if ~issparse(A) || ~issparse(B) || ~issparse(C)
                 error('LU factorization with full pivoting only works for sparse matrices.')
@@ -56,8 +56,8 @@
             obj.p = p;
             obj.q = q;
 
-            obj.k = k;
-            obj.t = t0+k;
+            obj.dt = dt;
+            obj.t = t0+dt;
             obj.n = 1;
             obj.v = v;
             obj.v_prev = v_prev;
@@ -69,7 +69,7 @@
         end
 
         function [vt,t] = getVt(obj)
-            vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
+            vt = (obj.v-obj.v_prev)/obj.dt; % Could be improved using u_tt = f(u))
             t = obj.t;
         end
 
@@ -77,7 +77,7 @@
         function E = getEnergy(obj)
             v  = obj.v;
             vp = obj.v_prev;
-            vt = (obj.v - obj.v_prev)/obj.k;
+            vt = (obj.v - obj.v_prev)/obj.dt;
 
             E = vt'*obj.A*vt + 1/2*(v'*obj.C*v + vp'*obj.C*vp);
         end
@@ -95,7 +95,7 @@
             obj.v(obj.q) = z;
 
             % Update time
-            obj.t = obj.t + obj.k;
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end