Mercurial > repos > public > sbplib

comparison +rv/+time/BDFDerivative.m @ 1014:e547794a9407 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add boot-strapping to RungeKuttaExteriorRV - Higher order BDF approximations are successively used as increasing number of time levels are obtained.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 06 Dec 2018 11:30:47 +0100
parents eb441fbdf379
children
comparison
equal deleted inserted replaced
1013:eb441fbdf379 1014:e547794a9407
1 class BDFDerivative < handle 1 classdef BDFDerivative < handle
2 properties 2 properties
3 dt 3 coefficients
4 coefficents
5 end 4 end
6 methods 5 methods %TBD: Decide if the BDF order should be passed in construction
7 function obj = BDFDerivative(order, dt) 6 % and only a row of coefficients stored based on the order.
8 obj.coefficents = obj.getBDFCoefficients(order); 7 % There would still be an implicit dependancy on the number
9 obj.dt = dt; 8 % of vectors in v_prev and elements in coefficients.
9 % In addition, dt could be stored, but this would be
10 % inflexible if different step sizes are employed.
11 function obj = BDFDerivative()
12 obj.coefficients = obj.getBDFCoefficients();
10 end 13 end
11 14 % Add asserts here?
12 function DvDt = evaluate(v, v_prev) 15 function DvDt = evaluate(obj, v, v_prev, dt)
13 DvDt = (obj.coefficients(1)*v - sum(obj.coefficients(2:end).*v_prev),2)/obj.dt; 16 order = size(v_prev,2);
17 DvDt = (obj.coefficients(order,1)*v - sum(obj.coefficients(order,2:order+1).*v_prev,2))/dt;
14 end 18 end
15 19 end
16 function c = getBDFCoefficients(obj, order) 20 methods(Static)
17 switch order 21 function c = getBDFCoefficients()
18 case 1 22 c = zeros(6,7);
19 c(1,1) = 1; c(1,2) = 1; 23 c(1,1) = 1; c(1,2) = 1;
20 case 2 24 c(2,1) = 3/2; c(2,2) = 4/2; c(2,3) = -1/2;
21 c(1,1) = 3/2; c(1,2) = 4/2; c(1,3) = -1/2; 25 c(3,1) = 11/6; c(3,2) = 18/6; c(3,3) = -9/6; c(3,4) = 2/6;
22 case 3 26 c(4,1) = 25/12; c(4,2) = 48/12; c(4,3) = -36/12; c(4,4) = 16/12; c(4,5) = -3/12;
23 c(1,1) = 11/6; c(1,2) = 18/6; c(1,3) = -9/6; c(1,4) = 2/6; 27 c(5,1) = 137/60; c(5,2) = 300/60; c(5,3) = -300/60; c(5,4) = 200/60; c(5,5) = -75/60; c(5,6) = 12/60;
24 case 4 28 c(6,1) = 147/60; c(6,2) = 360/60; c(6,3) = -450/60; c(6,4) = 400/60; c(6,5) = -225/60; c(6,6) = 72/60; c(6,7) = -10/60;
25 c(1,1) = 25/12; c(1,2) = 48/12; c(1,3) = -36/12; c(1,4) = 16/12; c(1,5) = -3/12;
26 case 5
27 c(1,1) = 137/60; c(1,2) = 300/60; c(1,3) = -300/60; c(1,4) = 200/60; c(1,5) = -75/60; c(1,6) = 12/60;
28 case 6
29 c(1,1) = 147/60; c(1,2) = 360/60; c(1,3) = -450/60; c(1,4) = 400/60; c(1,5) = -225/60; c(1,6) = 72/60; c(1,7) = -10/60;
30 assert(length(c) == (order+1));
31 end 29 end
32 end 30 end
33 end 31 end