Mercurial > repos > public > sbplib
diff +time/SBPInTimeSecondOrderFormImplicit.m @ 460:e0caae9ef6ed feature/grids
Add SBPinTime for linear DAE formulations (BUGS!)
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Mon, 19 Jun 2017 16:50:13 +0200 |
parents | +time/SBPInTimeSecondOrderForm.m@c92d2f8319c2 |
children | 45a3924140f4 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/SBPInTimeSecondOrderFormImplicit.m Mon Jun 19 16:50:13 2017 +0200 @@ -0,0 +1,80 @@ +classdef SBPInTimeSecondOrderFormImplicit < time.Timestepper + properties + A, B, C, f + AA, BB, ff + + n + t + k + + firstOrderTimeStepper + end + + methods + % Solves A*u_tt + B*u_t + C*u + f(t) = 0 + % A, B can either both be constants or both be function handles, + % They can also be omitted by setting them equal to the empty matrix. + function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, TYPE, order, blockSize) + default_arg('f', []); + default_arg('TYPE', []); + default_arg('order', []); + default_arg('blockSize',[]); + + m = length(v0); + + default_arg('A', sparse(m, m)); + default_arg('B', sparse(m, m)); + default_arg('C', sparse(m, m)); + + I = speye(m); + O = sparse(m,m); + + % Rewrite to + % AA*w_t = BB*w + ff(t); + + obj.AA = [ + I, O; + O, A; + ]; + obj.BB = [ + O, I; + -B, -C; + ]; + + if ~isempty(f) + obj.ff = @(t)[ + sparse(m,1); + -f(t); + ]; + else + obj.ff = @(t) sparse(2*m,1); + end + + w0 = [v0; v0t]; + + obj.k = k; + obj.t = t0; + obj.n = 0; + + obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize); + end + + function [v,t] = getV(obj) + w = obj.firstOrderTimeStepper.getV(); + v = w(1:end/2); + t = obj.t; + end + + function [vt,t] = getVt(obj) + w = obj.firstOrderTimeStepper.getV(); + vt = w(end/2+1:end); + t = obj.t; + end + + function obj = step(obj) + obj.firstOrderTimeStepper.step(); + obj.t = obj.firstOrderTimeStepper.t; + obj.n = obj.firstOrderTimeStepper.n; + end + end +end