Mercurial > repos > public > sbplib
changeset 992:bbd165cc585c feature/timesteppers
Move time.Rungekutta to time.rk.General
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 09 Jan 2019 10:59:38 +0100 |
| parents | a99f00896b8e |
| children | 44e7e497c3b7 |
| files | +time/+rk/General.m +time/Rungekutta.m |
| diffstat | 2 files changed, 81 insertions(+), 81 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/General.m Wed Jan 09 10:59:38 2019 +0100 @@ -0,0 +1,81 @@ +classdef General < time.Timestepper + properties + F % RHS of the ODE + dt % Time step + t % Time point + v % Solution vector + n % Time level + scheme % The scheme used for the time stepping, e.g rk4, rk6 etc. + coeffs % Butcher tableau coefficients + V % All stage approximations in most recent time step + K % All stage rates in most recent time step + end + + + methods + % Timesteps v_t = F(v,t), using the specified RK method from t = t0 with + % timestep dt and initial conditions v = v0 + function obj = General(F, dt, t0, v0, method, discreteData) + default_arg('method',"rk4"); + default_arg('discreteData', []); + obj.F = F; + obj.dt = dt; + obj.t = t0; + obj.v = v0; + obj.n = 0; + + % Extract the coefficients for the specified method + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = time.rk.butcherTableauFromStr(method); + obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); + + if isempty(discreteData) + 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 + end + + % v: Current solution + % t: Current time + % V: All stage approximations in most recent time step + % K: All stage rates in most recent time step + % T: Time points (corresponding to V and K) in most recent time step + function [v,t,V,T,K] = getV(obj) + v = obj.v; + t = obj.t; + V = obj.V; + K = obj.K; + T = obj.t + obj.dt*obj.coeffs.b; + end + + function obj = step(obj) + [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.n); + obj.t = obj.t + obj.dt; + obj.n = obj.n + 1; + end + + % Returns a vector of time points, including substage points, + % in the time interval [t0, tEnd]. + % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already. + function tvec = timePoints(obj, t0, tEnd) + N = round( (tEnd-t0)/obj.dt ); + tvec = zeros(N*obj.s, 1); + s = obj.coeffs.s; + c = obj.coeffs.c; + for i = 1:N + ind = (i-1)*s+1 : i*s; + tvec(ind) = ((i-1) + c')*obj.dt; + end + end + + % Returns a vector of quadrature weights corresponding to grid points + % in time interval [t0, tEnd], substage points included. + % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already. + function weights = quadWeights(obj, t0, tEnd) + N = round( (tEnd-t0)/obj.dt ); + b = obj.coeffs.b; + weights = repmat(b', N, 1); + end + end +end \ No newline at end of file
--- a/+time/Rungekutta.m Wed Jan 09 10:17:00 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,81 +0,0 @@ -classdef Rungekutta < time.Timestepper - properties - F % RHS of the ODE - dt % Time step - t % Time point - v % Solution vector - n % Time level - scheme % The scheme used for the time stepping, e.g rk4, rk6 etc. - coeffs % Butcher tableau coefficients - V % All stage approximations in most recent time step - K % All stage rates in most recent time step - end - - - methods - % Timesteps v_t = F(v,t), using the specified RK method from t = t0 with - % timestep dt and initial conditions v = v0 - function obj = Rungekutta(F, dt, t0, v0, method, discreteData) - default_arg('method',"rk4"); - default_arg('discreteData', []); - obj.F = F; - obj.dt = dt; - obj.t = t0; - obj.v = v0; - obj.n = 0; - - % Extract the coefficients for the specified method - % used for the RK updates from the Butcher tableua. - [s,a,b,c] = time.rk.butcherTableauFromStr(method); - obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); - - if isempty(discreteData) - 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 - end - - % v: Current solution - % t: Current time - % V: All stage approximations in most recent time step - % K: All stage rates in most recent time step - % T: Time points (corresponding to V and K) in most recent time step - function [v,t,V,T,K] = getV(obj) - v = obj.v; - t = obj.t; - V = obj.V; - K = obj.K; - T = obj.t + obj.dt*obj.coeffs.b; - end - - function obj = step(obj) - [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.n); - obj.t = obj.t + obj.dt; - obj.n = obj.n + 1; - end - - % Returns a vector of time points, including substage points, - % in the time interval [t0, tEnd]. - % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already. - function tvec = timePoints(obj, t0, tEnd) - N = round( (tEnd-t0)/obj.dt ); - tvec = zeros(N*obj.s, 1); - s = obj.coeffs.s; - c = obj.coeffs.c; - for i = 1:N - ind = (i-1)*s+1 : i*s; - tvec(ind) = ((i-1) + c')*obj.dt; - end - end - - % Returns a vector of quadrature weights corresponding to grid points - % in time interval [t0, tEnd], substage points included. - % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already. - function weights = quadWeights(obj, t0, tEnd) - N = round( (tEnd-t0)/obj.dt ); - b = obj.coeffs.b; - weights = repmat(b', N, 1); - end - end -end \ No newline at end of file