Mercurial > repos > public > sbplib
changeset 1100:27aaf8646a80 feature/timesteppers
Merge with default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 09 Apr 2019 21:48:30 +0200 |
| parents | d4fe089b2c4a (diff) 78d7e4e28e3e (current diff) |
| children | b895037bb701 |
| files | +scheme/Beam2d.m +scheme/TODO.txt +scheme/Wave2dCurve.m +scheme/error1d.m +scheme/error2d.m +scheme/errorMax.m +scheme/errorRelative.m +scheme/errorSbp.m +scheme/errorVector.m |
| diffstat | 12 files changed, 342 insertions(+), 354 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/ButcherTableau.m Tue Apr 09 21:48:30 2019 +0200 @@ -0,0 +1,127 @@ +classdef ButcherTableau + properties + a,b,c + end + + methods + % A ButcherTableau describes a specific rungekutta method where + % y(n+1) = y(n) + dt*b(i)*k(i) + % k(i) = F(t + c(i)*dt, Y(i)) + % Y(i) = y(i) + dt*a(i,j)*k(j) + % where repeating indecies imply summation + function obj = ButcherTableau(a,b,c) + s = length(c); + assertSize(a, [s,s]); + assertLength(b, s); + + obj.a = a; + obj.b = b; + obj.c = c; + end + + function s = nStages(obj) + s = length(obj.c); + end + + function b = isExplicit(obj) + b = all(all(triu(obj.a)==0)); + end + + function g = testEquationGain(obj, z) + default_arg('z', sym('z')); + s = obj.nStages(); + + b = sym(obj.b); + A = sym(obj.a); + one = sym(ones(s,1)); + I = sym(eye(s)); + + g = abs(1 + z*b*inv(I-z*A)*one); + end + + % TBD: Add functions for checking accuracy, stability? + end + + methods(Static) + % TVD (Total Variational Diminishing) + function bt = tvd_3() + a = [ + 0, 0, 0; + 1, 0, 0; + 1/4, 1/4, 0; + ]; + b = [1/6, 1/6, 2/3]; + c = [0 1 1/2]; + + bt = time.rk.ButcherTableau(a,b,c); + end + + % Standard RK4 + function bt = rk4() + a = [ + 0, 0, 0, 0; + 1/2, 0, 0, 0; + 0, 1/2, 0, 0; + 0, 0, 1, 0; + ]; + + b = [1/6 1/3 1/3 1/6]; + c = [0, 1/2, 1/2, 1]; + + bt = time.rk.ButcherTableau(a,b,c); + end + + % 3/8 RK4 (Kuttas method). Lower truncation error, more flops. + % Irreducible, unlike standard rk4. + function bt = rk4_3_8() + a = [ + 0, 0, 0, 0; + 1/3, 0, 0, 0; + -1/3, 1, 0, 0; + 1, -1, 1, 0; + ]; + + b = [1/8 3/8 3/8 1/8]; + c = [0, 1/3, 2/3, 1]; + + bt = time.rk.ButcherTableau(a,b,c); + end + + % Runge-Kutta 6 from Alshina07 + function bt = rk6() + a = zeros(7,7); + + a(2,1) = 4/7; + + a(3,1) = 115/112; + a(3,2) = -5/16; + + a(4,1) = 589/630; + a(4,2) = 5/18; + a(4,3) = -16/45; + + a(5,1) = 229/1200 - 29/6000*sqrt(5); + a(5,2) = 119/240 - 187/1200*sqrt(5); + a(5,3) = -14/75 + 34/375*sqrt(5); + a(5,4) = -3/100*sqrt(5); + + a(6,1) = 71/2400 - 587/12000*sqrt(5); + a(6,2) = 187/480 - 391/2400*sqrt(5); + a(6,3) = -38/75 + 26/375*sqrt(5); + a(6,4) = 27/80 - 3/400*sqrt(5); + a(6,5) = (1+sqrt(5))/4 + + a(7,1) = -49/480 + 43/160*sqrt(5); + a(7,2) = -425/96 + 51/32*sqrt(5); + a(7,3) = 52/15 - 4/5*sqrt(5); + a(7,4) = -27/16 + 3/16*sqrt(5); + a(7,5) = 5/4 - 3/4*sqrt(5); + a(7,6) = 5/2 - 1/2*sqrt(5); + + b = [1/12 0 0 0 5/12 5/12 1/12]; + c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1]; + + bt = time.rk.ButcherTableau(a,b,c); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/Explicit.m Tue Apr 09 21:48:30 2019 +0200 @@ -0,0 +1,142 @@ +classdef Explicit < time.Timestepper + properties + F % RHS of the ODE + dt % Time step + t % Time point + v % Solution vector + n % Time level + bt + end + + + methods + % Timesteps v_t = F(t,v), using the specified ButcherTableau + % from t = t0 with timestep dt and initial conditions v(0) = v0 + function obj = Explicit(F, dt, t0, v0, bt) + assertType(bt, 'time.rk.ButcherTableau') + obj.F = F; + obj.dt = dt; + obj.t = t0; + obj.v = v0; + obj.n = 0; + + assert(bt.isExplicit()) + obj.bt = bt; + 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] = getV(obj) + v = obj.v; + t = obj.t; + end + + function obj = step(obj) + s = obj.bt.nStages(); + a = obj.bt.a; + b = obj.bt.b; + c = obj.bt.c; + + % Compute rates K + K = zeros(length(v), s); + for i = 1:s + V_i = obj.v; + for j = 1:i-1 + V_i = V_i + dt*a(i,j)*K(:,j); + end + K(:,i) = F(t+dt*c(i), V_i); + end + + % Compute updated solution + v_next = v; + for i = 1:s + v_next = v_next + dt*b(i)*K(:,i); + end + + obj.v = v_next; + obj.t = obj.t + obj.dt; + obj.n = obj.n + 1; + end + + % TBD: Method name + % TBD: Parameter name + % + % Takes a regular step but with discreteRates(:,i) added to RHS for stage i. + % v_t = F(t,v) + discreteRates(:, ...) + % + % Also returns the stage approximations (V) and stage rates (K). + function [v,t, V, K] = stepWithDiscreteData(obj, discreteRates) + s = obj.bt.nStages(); + a = obj.bt.a; + b = obj.bt.b; + c = obj.bt.c; + + % Compute rates K and stage approximations V + K = zeros(length(v), s); + V = zeros(length(v), s); + for i = 1:s + V_i = obj.v; + for j = 1:i-1 + V_i = V_i + dt*a(i,j)*K(:,j); + end + + K_i = F(t+dt*c(i), V_i); + K_i = K_i + discreteRates(:,i); + + V(:,i) = V_i; + K(:,i) = K_i; + end + + % Compute updated updated solution + v_next = v; + for i = 1:s + v_next = v_next + dt*b(i)*K(:,i); + end + + obj.v = v_next; + 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) + % TBD: Should this be implemented here or somewhere else? + 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) + % TBD: Should this be implemented here or somewhere else? + N = round( (tEnd-t0)/obj.dt ); + b = obj.coeffs.b; + weights = repmat(b', N, 1); + end + end + + methods(Static) + % TBD: Function name + function ts = methodFromStr(F, dt, t0, v0, methodStr) + try + bt = time.rk.ButcherTableau.(method); + catch + error('Runge-Kutta method ''%s'' is not implemented', methodStr) + end + + ts = time.rk.Explicit(F, dt, t0, v0, bt); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/get_rk4_time_step.m Tue Apr 09 21:48:30 2019 +0200 @@ -0,0 +1,21 @@ +% Calculate the size of the largest time step given the largest evalue for a operator with pure imaginary e.values. +function k = get_rk4_time_step(lambda,l_type) + default_arg('l_type','complex') + + rad = abs(lambda); + if strcmp(l_type,'real') + % Real eigenvalue + % kl > -2.7852 + k = 2.7852/rad; + + elseif strcmp(l_type,'imag') + % Imaginary eigenvalue + % |kl| < 2.8284 + k = 2.8284/rad; + elseif strcmp(l_type,'complex') + % |kl| < 2.5 + k = 2.5/rad; + else + error('l_type must be one of ''real'',''imag'' or ''complex''.') + end +end \ No newline at end of file
--- a/+time/+rk4/get_rk4_time_step.m Mon Apr 08 20:15:37 2019 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,21 +0,0 @@ -% Calculate the size of the largest time step given the largest evalue for a operator with pure imaginary e.values. -function k = get_rk4_time_step(lambda,l_type) - default_arg('l_type','complex') - - rad = abs(lambda); - if strcmp(l_type,'real') - % Real eigenvalue - % kl > -2.7852 - k = 2.7852/rad; - - elseif strcmp(l_type,'imag') - % Imaginary eigenvalue - % |kl| < 2.8284 - k = 2.8284/rad; - elseif strcmp(l_type,'complex') - % |kl| < 2.5 - k = 2.5/rad; - else - error('l_type must be one of ''real'',''imag'' or ''complex''.') - end -end \ No newline at end of file
--- a/+time/+rk4/rk4_stability.m Mon Apr 08 20:15:37 2019 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,58 +0,0 @@ -function rk_stability() - ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4)); - circ = @(z)(abs(z)); - - - % contour(X,Y,z) - ax = [-4 2 -3 3]; - % hold on - fcontour(ruku4,[1,1],[-3, 0.6],[-3.2, 3.2]) - hold on - r = 2.6; - fcontour(circ,[r,r],[-3, 0.6],[-3.2, 3.2],'r') - hold off - % contour(X,Y,z,[1,1],'b') - axis(ax) - title('4th order Runge-Kutta stability region') - xlabel('Re') - ylabel('Im') - axis equal - grid on - box on - hold off - % surf(X,Y,z) - - - rk4roots() -end - -function fcontour(f,levels,x_lim,y_lim,opt) - default_arg('opt','b') - x = linspace(x_lim(1),x_lim(2)); - y = linspace(y_lim(1),y_lim(2)); - [X,Y] = meshgrid(x,y); - mu = X+ 1i*Y; - - z = f(mu); - - contour(X,Y,z,levels,opt) - -end - - -function rk4roots() - ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4)); - % Roots for real evalues: - F = @(x)(abs(ruku4(x))-1); - real_x = fzero(F,-3); - - % Roots for imaginary evalues: - F = @(x)(abs(ruku4(1i*x))-1); - imag_x1 = fzero(F,-3); - imag_x2 = fzero(F,3); - - - fprintf('Real x = %f\n',real_x) - fprintf('Imag x = %f\n',imag_x1) - fprintf('Imag x = %f\n',imag_x2) -end \ No newline at end of file
--- a/+time/+rk4/rungekutta_4.m Mon Apr 08 20:15:37 2019 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -% Takes one time step of size k using the rungekutta method -% starting from v_0 and where the function F(v,t) gives the -% time derivatives. -function v = rungekutta_4(v, t , k, F) - k1 = F(v ,t ); - k2 = F(v+0.5*k*k1,t+0.5*k); - k3 = F(v+0.5*k*k2,t+0.5*k); - k4 = F(v+ k*k3,t+ k); - v = v + (1/6)*(k1+2*(k2+k3)+k4)*k; -end \ No newline at end of file
--- a/+time/+rk4/rungekutta_6.m Mon Apr 08 20:15:37 2019 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% Takes one time step of size k using the rungekutta method -% starting from v_0 and where the function F(v,t) gives the -% time derivatives. -function v = rungekutta_6(v, t , k, F) - s = 7 - k = zeros(length(v),s) - a = zeros(7,6); - c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1]; - b = [1/12, 0, 0, 0, 5/12, 5/12, 1/12]; - a = [ - 0, 0, 0, 0, 0, 0; - 4/7, 0, 0, 0, 0, 0; - 115/112, -5/16, 0, 0, 0, 0; - 589/630, 5/18, -16/45, 0, 0, 0; - 229/1200 - 29/6000*sqrt(5), 119/240 - 187/1200*sqrt(5), -14/75 + 34/375*sqrt(5), -3/100*sqrt(5), 0, 0; - 71/2400 - 587/12000*sqrt(5), 187/480 - 391/2400*sqrt(5), -38/75 + 26/375*sqrt(5), 27/80 - 3/400*sqrt(5), (1+sqrt(5))/4, 0; - -49/480 + 43/160*sqrt(5), -425/96 + 51/32*sqrt(5), 52/15 - 4/5*sqrt(5), -27/16 + 3/16*sqrt(5), 5/4 - 3/4*sqrt(5), 5/2 - 1/2*sqrt(5); - ] - - for i = 1:s - u = v - for j = 1: i-1 - u = u + h*a(i,j) * k(:,j) - end - k(:,i) = F(t+c(i)*k,u) - end - - for i = 1:s - v = v + k*b(i)*k(:,i) - end -end
--- a/+time/Rk4SecondOrderNonlin.m Mon Apr 08 20:15:37 2019 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,77 +0,0 @@ -classdef Rk4SecondOrderNonlin < time.Timestepper - properties - F - k - t - w - m - - D - E - S - - n - end - - - methods - function obj = Rk4SecondOrderNonlin(D, E, S, k, t0, v0, v0t) - default_arg('S',0); - default_arg('E',0); - - if isnumeric(S) - S = @(v,t)S; - end - - if isnumeric(E) - E = @(v)E; - end - - obj.k = k; - obj.t = t0; - obj.w = [v0; v0t]; - - m = length(v0); - function wt = F(w,t) - v = w(1:m); - vt = w(m+1:end); - - % Def: w = [v; vt] - wt(1:m,1) = vt; - wt(m+1:2*m,1) = D(v)*v + E(v)*vt + S(v,t); - - end - - obj.F = @F; - obj.D = D; - obj.E = E; - obj.S = S; - obj.m = m; - obj.n = 0; - end - - function [v,t] = getV(obj) - v = obj.w(1:end/2); - t = obj.t; - end - - function [vt,t] = getVt(obj) - vt = obj.w(end/2+1:end); - t = obj.t; - end - - function obj = step(obj) - obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); - obj.t = obj.t + obj.k; - obj.n = obj.n + 1; - end - end - - - methods (Static) - function k = getTimeStep(lambda) - k = rk4.get_rk4_time_step(lambda); - end - end - -end \ No newline at end of file
--- a/+time/Rungekutta4.m Mon Apr 08 20:15:37 2019 +0000 +++ b/+time/Rungekutta4.m Tue Apr 09 21:48:30 2019 +0200 @@ -1,36 +1,23 @@ classdef Rungekutta4 < time.Timestepper properties - D - S F - k + dt t v - m n end methods - function obj = Rungekutta4(D, S, k, t0, v0) - obj.D = D; - obj.k = k; + % Create a time stepper for + % v_t = F(t,v), v(t0) = v0 + % with step size dt. + function obj = Rungekutta4(F, dt, t0, v0) + obj.F = F; + obj.dt = dt; obj.t = t0; obj.v = v0; - obj.m = length(v0); obj.n = 0; - - if S == 0 - obj.S = zeros(obj.m,1); - else - obj.S = S; - end - - if S == 0 - obj.F = @(v,t)(obj.D*v); - else - obj.F = @(v,t)(obj.D*v + obj.S); - end end function [v,t] = getV(obj) @@ -39,17 +26,24 @@ end function obj = step(obj) - obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F); - obj.t = obj.t + obj.k; + t = obj.t; + v = obj.v; + dt = obj.dt; + + k1 = obj.F(t, v); + k2 = obj.F(t + 0.5*dt, v + 0.5*dt*k1); + k3 = obj.F(t + 0.5*dt, v + 0.5*dt*k2); + k4 = obj.F(t + dt, v + dt*k3); + + obj.v = v + dt*(1/6)*(k1+2*(k2+k3)+k4); + obj.t = obj.t + obj.dt; obj.n = obj.n + 1; end end - methods (Static) - function k = getTimeStep(lambda) - k = rk4.get_rk4_time_step(lambda); + function dt = getTimeStep(lambda) + dt = rk4.get_rk4_time_step(lambda); end end - end \ No newline at end of file
--- a/+time/Rungekutta4SecondOrder.m Mon Apr 08 20:15:37 2019 +0000 +++ b/+time/Rungekutta4SecondOrder.m Tue Apr 09 21:48:30 2019 +0200 @@ -1,115 +1,57 @@ classdef Rungekutta4SecondOrder < time.Timestepper properties F - k - t - w - m - D - E - S - M - C - n + dt + t, n + v, v_t end methods - % Solves u_tt = Du + Eu_t + S by - % Rewriting on first order form: - % w_t = M*w + C(t) - % where - % M = [ - % 0, I; - % D, E; - % ] - % and - % C(t) = [ - % 0; - % S(t) - % ] - % D, E, S can either all be constants or all be function handles, - % They can also be omitted by setting them equal to the empty matrix. - function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) - obj.D = D; - obj.E = E; - obj.S = S; - obj.m = length(v0); + % Create a time stepper for + % v_tt = F(t,v,v_t), v(t0) = v0, v_t(t0) = v0t + % with step size dt, by rewriting on first order form + function obj = Rungekutta4SecondOrder(F, dt, t0, v0, v0t) + obj.F = F; + obj.dt = dt; + obj.t = t0; obj.n = 0; - - if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle') - default_arg('D', @(t)sparse(obj.m, obj.m)); - default_arg('E', @(t)sparse(obj.m, obj.m)); - default_arg('S', @(t)sparse(obj.m, 1) ); - - if ~isa(D, 'function_handle') - D = @(t)D; - end - if ~isa(E, 'function_handle') - E = @(t)E; - end - if ~isa(S, 'function_handle') - S = @(t)S; - end - - obj.k = k; - obj.t = t0; - obj.w = [v0; v0t]; - - % Avoid matrix formulation because it is VERY slow - obj.F = @(w,t)[ - w(obj.m+1:end); - D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t); - ]; - else - - default_arg('D', sparse(obj.m, obj.m)); - default_arg('E', sparse(obj.m, obj.m)); - default_arg('S', sparse(obj.m, 1) ); - - I = speye(obj.m); - O = sparse(obj.m,obj.m); - - obj.M = [ - O, I; - D, E; - ]; - obj.C = [ - zeros(obj.m,1); - S; - ]; - - obj.k = k; - obj.t = t0; - obj.w = [v0; v0t]; - - obj.F = @(w,t)(obj.M*w + obj.C); - end + obj.v = v0; + obj.v_t = v0t; end function [v,t] = getV(obj) - v = obj.w(1:end/2); + v = obj.v t = obj.t; end function [vt,t] = getVt(obj) - vt = obj.w(end/2+1:end); + vt = obj.v_t; t = obj.t; end function obj = step(obj) - obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); - obj.t = obj.t + obj.k; + t = obj.t; + v = obj.v; + v_t = obj.v_t; + dt = obj.dt; + + k1 = obj.F(t, v, v_t); + k2 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t, v_t + 1/2*dt*k1); + k3 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t + 1/4*dt^2*k1, v_t + 1/2*dt*k2); + k4 = obj.F(t + dt, v + dt*v_t + 1/2*dt^2*k2, v_t + dt*k3); + + obj.v = v + dt*v_t + dt^2*(1/6)*(k1 + k2 + k3); + obj.v_t = v_t + dt*(1/6)*(k1 + 2*k2 + 2*k3 + k4); + obj.t = obj.t + obj.dt; obj.n = obj.n + 1; end end - methods (Static) function k = getTimeStep(lambda) k = rk4.get_rk4_time_step(lambda); end end - -end \ No newline at end of file +end
--- a/+time/Rungekutta4proper.m Mon Apr 08 20:15:37 2019 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,42 +0,0 @@ -classdef Rungekutta4proper < time.Timestepper - properties - F - k - t - v - m - n - end - - - methods - % Timesteps v_t = F(v,t), using RK4 fromt t = t0 with timestep k and initial conditions v = v0 - function obj = Rungekutta4proper(F, k, t0, v0) - obj.F = F; - obj.k = k; - obj.t = t0; - obj.v = v0; - obj.m = length(v0); - obj.n = 0; - end - - function [v,t] = getV(obj) - v = obj.v; - t = obj.t; - end - - function obj = step(obj) - obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F); - obj.t = obj.t + obj.k; - obj.n = obj.n + 1; - end - end - - - methods (Static) - function k = getTimeStep(lambda) - k = rk4.get_rk4_time_step(lambda); - end - end - -end \ No newline at end of file
--- a/+time/Timestepper.m Mon Apr 08 20:15:37 2019 +0000 +++ b/+time/Timestepper.m Tue Apr 09 21:48:30 2019 +0200 @@ -1,13 +1,14 @@ classdef Timestepper < handle properties (Abstract) t - k + k % TBD: should this be a method instead? n end methods (Abstract) - [v,t] = getV(obj) - obj = step(obj) + % Returns the solution vector v at timestep t. + [v,t] = getV(obj) + obj = step(obj) end