Mercurial > repos > public > sbplib
changeset 888:8732d6bd9890 feature/timesteppers
Add general Runge-Kutta class
- Add a general Runge-Kutta class which time integrates the solution based on coefficients obtained from a Butcher tableau
- Add butcher tableau which returns coefficents for the specified Runge-Kutta method
- Remove RungKutta4proper, since obsolete
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 15 Nov 2018 17:10:01 -0800 |
| parents | 50d5a3843099 |
| children | f5e14e5986b5 |
| files | +time/+rk/butcherTableau.m +time/+rk/rungekutta.m +time/+rk/rungekutta_6.m +time/Rungekutta.m +time/Rungekutta4proper.m |
| diffstat | 5 files changed, 116 insertions(+), 78 deletions(-) [+] |
line wrap: on
line diff
diff -r 50d5a3843099 -r 8732d6bd9890 +time/+rk/butcherTableau.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/butcherTableau.m Thu Nov 15 17:10:01 2018 -0800 @@ -0,0 +1,50 @@ +% Returns the coefficients used in a RK method as defined by a Butcher Tableau. +% +% @param method - string specifying which Runge-Kutta method to be used. +% @return s - number of stages +% @return a - coefficents for intermediate stages +% @return b - weights for summing stages +% @return c - time step coefficents for intermediate stages +function [s,a,b,c] = butcherTableau(method) +switch method + case "tvd-3" + % TVD (Total Variational Diminishing) + s = 3; + a = zeros(s,s-1); + a(2,1) = 1; + a(3,1) = 1/4; a(3,2) = 1/4; + b = [1/6, 1/6, 2/3]; + c = [0 1 1/2]; + case "rk4" + % Standard RK4 + s = 4; + a = zeros(s,s-1); + a(2,1) = 1/2; + a(3,1) = 0; a(3,2) = 1/2; + a(4,1) = 0; a(4,2) = 0; a(4,3) = 1; + b = [1/6 1/3 1/3 1/6]; + c = [0, 1/2, 1/2, 1]; + case "rk4-3/8" + % 3/8 RK4 (Kuttas method). Lower truncation error, more flops + s = 4; + a = zeros(s,s-1); + a(2,1) = 1/3; + a(3,1) = -1/3; a(3,2) = 1; + a(4,1) = 1; a(4,2) = -1; a(4,3) = 1; + b = [1/8 3/8 3/8 1/8]; + c = [0, 1/3, 2/3, 1]; + case "rk6" + % Runge-Kutta 6 from Alshina07 + s = 7; + a = zeros(s,s-1); + 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]; + otherwise + error('That Runge-Kutta method is not implemented', method) +end \ No newline at end of file
diff -r 50d5a3843099 -r 8732d6bd9890 +time/+rk/rungekutta.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rungekutta.m Thu Nov 15 17:10:01 2018 -0800 @@ -0,0 +1,20 @@ +% Takes one time step of size dt using the rungekutta method +% starting from @arg v and where the function F(v,t) gives the +% time derivatives. coeffs is a struct holding the RK coefficients +% for the specific method. +function v = rungekutta(v, t , dt, F, coeffs) + % Compute the intermediate stages k + k = zeros(length(v), coeffs.s); + for i = 1:coeffs.s + u = v; + for j = 1:i-1 + u = u + dt*coeffs.a(i,j)*k(:,j); + end + k(:,i) = F(u,t+coeffs.c(i)*dt); + end + % Compute the updated solution as a linear combination + % of the intermediate stages. + for i = 1:coeffs.s + v = v + dt*coeffs.b(i)*k(:,i); + end +end \ No newline at end of file
diff -r 50d5a3843099 -r 8732d6bd9890 +time/+rk/rungekutta_6.m --- a/+time/+rk/rungekutta_6.m Thu Nov 15 16:42:58 2018 -0800 +++ /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
diff -r 50d5a3843099 -r 8732d6bd9890 +time/Rungekutta.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Rungekutta.m Thu Nov 15 17:10:01 2018 -0800 @@ -0,0 +1,46 @@ +classdef Rungekutta < time.Timestepper + properties + F % RHS of the ODE + k % Time step + t % Time point + v % Solution vector + n % Time level + scheme % The scheme used for the time stepping, e.g rk4, rk6 etc. + end + + + methods + % Timesteps v_t = F(v,t), using the specified RK method from t = t0 with + % timestep k and initial conditions v = v0 + function obj = Rungekutta(F, k, t0, v0, method) + default_arg('method',"rk4"); + obj.F = F; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + % TODO: method "rk4" is also implemented in the butcher tableau, but the rungekutta_4.m implementation + % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it. + if (method == "rk4") + obj.scheme = @time.rk.rungekutta_4; + else + % Extract the coefficients for the specified method + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = time.rk.butcherTableau(method); + coeffs = struct('s',s,'a',a,'b',b,'c',c); + obj.scheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs); + end + end + + function [v,t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function obj = step(obj) + obj.v = obj.scheme(obj.v, obj.t, obj.k, obj.F); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file
diff -r 50d5a3843099 -r 8732d6bd9890 +time/Rungekutta4proper.m --- a/+time/Rungekutta4proper.m Thu Nov 15 16:42:58 2018 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +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 state = getState(obj) - state = struct('v', obj.v, 't', obj.t, 'k', obj.k); - end - - function obj = step(obj) - obj.v = time.rk.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