Mercurial > repos > public > sbplib

comparison +time/+rk/General.m @ 992:bbd165cc585c feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Move time.Rungekutta to time.rk.General
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 09 Jan 2019 10:59:38 +0100
parents +time/Rungekutta.m@a32856fc2ad2
children 44e7e497c3b7
comparison
equal deleted inserted replaced
991:a99f00896b8e 992:bbd165cc585c
1 classdef General < time.Timestepper
2 properties
3 F % RHS of the ODE
4 dt % Time step
5 t % Time point
6 v % Solution vector
7 n % Time level
8 scheme % The scheme used for the time stepping, e.g rk4, rk6 etc.
9 coeffs % Butcher tableau coefficients
10 V % All stage approximations in most recent time step
11 K % All stage rates in most recent time step
12 end
13
14
15 methods
16 % Timesteps v_t = F(v,t), using the specified RK method from t = t0 with
17 % timestep dt and initial conditions v = v0
18 function obj = General(F, dt, t0, v0, method, discreteData)
19 default_arg('method',"rk4");
20 default_arg('discreteData', []);
21 obj.F = F;
22 obj.dt = dt;
23 obj.t = t0;
24 obj.v = v0;
25 obj.n = 0;
26
27 % Extract the coefficients for the specified method
28 % used for the RK updates from the Butcher tableua.
29 [s,a,b,c] = time.rk.butcherTableauFromStr(method);
30 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
31
32 if isempty(discreteData)
33 obj.scheme = @(v,t,n) time.rk.rungekutta(v, t, dt, F, obj.coeffs);
34 else
35 obj.scheme = @(v,t,n) time.rk.rungekuttaDiscreteData(v, t, dt, F, obj.coeffs, discreteData, n);
36 end
37 end
38
39 % v: Current solution
40 % t: Current time
41 % V: All stage approximations in most recent time step
42 % K: All stage rates in most recent time step
43 % T: Time points (corresponding to V and K) in most recent time step
44 function [v,t,V,T,K] = getV(obj)
45 v = obj.v;
46 t = obj.t;
47 V = obj.V;
48 K = obj.K;
49 T = obj.t + obj.dt*obj.coeffs.b;
50 end
51
52 function obj = step(obj)
53 [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.n);
54 obj.t = obj.t + obj.dt;
55 obj.n = obj.n + 1;
56 end
57
58 % Returns a vector of time points, including substage points,
59 % in the time interval [t0, tEnd].
60 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
61 function tvec = timePoints(obj, t0, tEnd)
62 N = round( (tEnd-t0)/obj.dt );
63 tvec = zeros(N*obj.s, 1);
64 s = obj.coeffs.s;
65 c = obj.coeffs.c;
66 for i = 1:N
67 ind = (i-1)*s+1 : i*s;
68 tvec(ind) = ((i-1) + c')*obj.dt;
69 end
70 end
71
72 % Returns a vector of quadrature weights corresponding to grid points
73 % in time interval [t0, tEnd], substage points included.
74 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
75 function weights = quadWeights(obj, t0, tEnd)
76 N = round( (tEnd-t0)/obj.dt );
77 b = obj.coeffs.b;
78 weights = repmat(b', N, 1);
79 end
80 end
81 end