Mercurial > repos > public > sbplib
changeset 846:c6fcee3fcf1b feature/burgers1d
Add generalized RungeKutta and RungeKuttaRV class which extracts its coefficients from a butcher tableau, specified on the scheme.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 20 Sep 2018 17:51:19 +0200 |
| parents | 1e057b0f2fed |
| children | 1c6f1595bb94 |
| files | +time/+rk/butcherTableau.m +time/+rk/get_rk4_time_step.m +time/+rk/rk4_stability.m +time/+rk/rungekutta.m +time/+rk/rungekuttaRV.m +time/+rk/rungekutta_4.m +time/+rk/rungekutta_4RV.m +time/+rk/rungekutta_6.m +time/+rk/rungekutta_6RV.m +time/+rk4/get_rk4_time_step.m +time/+rk4/rk4_stability.m +time/+rk4/rungekutta_4.m +time/+rk4/rungekutta_4RV.m +time/+rk4/rungekutta_6.m +time/+rk4/rungekutta_6RV.m +time/Rk4SecondOrderNonlin.m +time/Rungekutta.m +time/Rungekutta4.m +time/Rungekutta4RV.m +time/Rungekutta4SecondOrder.m +time/Rungekutta4proper.m +time/RungekuttaRV.m |
| diffstat | 22 files changed, 356 insertions(+), 242 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/butcherTableau.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,37 @@ +function [s,a,b,c] = butcherTableau(order) + +switch order + + case 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 4 + % 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 6 + % 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 order is not implemented', order) + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/get_rk4_time_step.m Thu Sep 20 17:51:19 2018 +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rk4_stability.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,58 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rungekutta.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,20 @@ +% Takes one time step of size dt using the rungekutta method +% starting from v_0 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rungekuttaRV.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,31 @@ +% Takes one time step of size dt using the rungekutta method +% starting from v_0 and where the function F(v,t,RV) gives the +% time derivatives. coeffs is a struct holding the RK coefficients +% for the specific method. RV is the residual viscosity which is updated +% in between the stages and after the updated solution is computed. +function v = rungekuttaRV(v, t , dt, F, RV, coeffs) + + % Move one stage outside to avoid branching for updating the + % residual inside the loop. + k = zeros(length(v), coeffs.s); + k(:,1) = F(v,t,RV.getViscosity()); + + % Compute the intermediate stages k + for i = 2:coeffs.s + u = v; + for j = 1:i-1 + u = u + dt*coeffs.a(i,j)*k(:,j); + end + RV.update(u,v,coeffs.c(i)*dt); + k(:,i) = F(u,t+coeffs.c(i)*dt, RV.getViscosity()); + end + + % Compute the updated solution as a linear combination + % of the intermediate stages. + u = v; + for i = 1:coeffs.s + u = u + dt*coeffs.b(i)*k(:,i); + end + RV.update(u,v,dt); + v = u; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rungekutta_4.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,10 @@ +% Takes one time step of size dt 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 , dt, F) + k1 = F(v ,t ); + k2 = F(v+0.5*dt*k1,t+0.5*dt); + k3 = F(v+0.5*dt*k2,t+0.5*dt); + k4 = F(v+ dt*k3,t+ dt); + v = v + (1/6)*(k1+2*(k2+k3)+k4)*dt; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rungekutta_4RV.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,18 @@ +% 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_4RV(v, t , k, F, RV) + k1 = F(v, t, RV.getViscosity()); + + RV.update(v+0.5*k*k1, v, 0.5*k); + k2 = F(v+0.5*k*k1, t+0.5*k, RV.getViscosity()); + + RV.update(v+0.5*k*k2, v, 0.5*k); + k3 = F(v+0.5*k*k2, t+0.5*k, RV.getViscosity()); + + RV.update(v+k*k3, v, k); + k4 = F(v+k*k3,t+k, RV.getViscosity()); + + RV.update(v + (1/6)*(k1+2*(k2+k3)+k4)*k, v, k); + v = v + (1/6)*(k1+2*(k2+k3)+k4)*k; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rungekutta_6.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,31 @@ +% Takes one time step of size dt 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 , dt, 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 + dt*a(i,j)*k(:,j); + end + k(:,i) = F(t+c(i)*dt,u); + end + + for i = 1:s + v = v + dt*b(i)*k(:,i); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rungekutta_6RV.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,37 @@ +% Takes one time step of size dt using the rungekutta method +% starting from v_0 and where the function F(v,t) gives the +% time derivatives. +function v = rungekutta_6RV(v, t , dt, F, RV) + 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); + ]; + + k(:,1) = F(v,t,RV.getViscosity()); + + for i = 2:s + u = v; + for j = 1:i-1 + u = u + dt*a(i,j)*k(:,j); + end + RV.update(u,v,c(i)*dt); + k(:,i) = F(u,t+c(i)*dt,RV.getViscosity()); + end + + u = v; + for i = 1:s + u = u + dt*b(i)*k(:,i); + end + RV.update(u,v,dt); + v = u; +end
--- a/+time/+rk4/get_rk4_time_step.m Wed Sep 19 16:32:05 2018 +0200 +++ /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 Wed Sep 19 16:32:05 2018 +0200 +++ /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 Wed Sep 19 16:32:05 2018 +0200 +++ /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_4RV.m Wed Sep 19 16:32:05 2018 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,34 +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_4RV(v, t , k, F, RV) - - v1 = v + k/2*F(v,t, RV.getViscosity()); - - RV.update(v1, v, k/2); - v2 = v + k/2*F(v1,t+k/2, RV.getViscosity()); - - RV.update(v2, v, k/2); - v3 = v + k*F(v2,t+k/2, RV.getViscosity()); - - RV.update(v3,v,k); - v4 = v + k*F(v3,t+k, RV.getViscosity()); - - v_next = 1/6*(-3*v + 2*v1 + 4*v2 + 2*v3 + v4); - RV.update(v_next,v,k); - v = v_next; - - % k1 = F(v, t, RV.getViscosity()); - - % RV.update(v+0.5*k*k1, v, 0.5*k); - % k2 = F(v+0.5*k*k1, t+0.5*k, RV.getViscosity()); - - % RV.update(v+0.5*k*k2, v, 0.5*k); - % k3 = F(v+0.5*k*k2, t+0.5*k, RV.getViscosity()); - - % RV.update(v+k*k3, v, k); - % k4 = F(v+k*k3,t+k, RV.getViscosity()); - - % RV.update(v + (1/6)*(k1+2*(k2+k3)+k4)*k, v, k); - % v = v + (1/6)*(k1+2*(k2+k3)+k4)*k; -end \ No newline at end of file
--- a/+time/+rk4/rungekutta_6.m Wed Sep 19 16:32:05 2018 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% Takes one time step of size dt 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 , dt, 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 + dt*a(i,j)*k(:,j); - end - k(:,i) = F(t+c(i)*dt,u); - end - - for i = 1:s - v = v + dt*b(i)*k(:,i); - end -end
--- a/+time/+rk4/rungekutta_6RV.m Wed Sep 19 16:32:05 2018 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,34 +0,0 @@ -% Takes one time step of size dt using the rungekutta method -% starting from v_0 and where the function F(v,t) gives the -% time derivatives. -function v = rungekutta_6RV(v, t , dt, F, RV) - 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 + dt*a(i,j)*k(:,j); - end - if (i > 1) - RV.update(u,v,c(i)*dt); - end - k(:,i) = F(u,t+c(i)*dt,RV.getViscosity()); - end - - for i = 1:s - v = v + dt*b(i)*k(:,i); - end -end
--- a/+time/Rk4SecondOrderNonlin.m Wed Sep 19 16:32:05 2018 +0200 +++ b/+time/Rk4SecondOrderNonlin.m Thu Sep 20 17:51:19 2018 +0200 @@ -61,7 +61,7 @@ end function obj = step(obj) - obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); + obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F); obj.t = obj.t + obj.k; obj.n = obj.n + 1; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Rungekutta.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,39 @@ +classdef Rungekutta < time.Timestepper + properties + F % RHS of the ODE + k % Time step + t % Time point + v % Solution vector + n % Time level + coeffs % The coefficents used for the RK time integration + end + + + methods + % Timesteps v_t = F(v,t), using RK with specfied order from t = t0 with + % timestep k and initial conditions v = v0 + function obj = Rungekutta(F, k, t0, v0, order) + default_arg('order','4'); + obj.F = F; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + % Extract the coefficients for the specified order + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = butcherTableau(order); + obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); + end + + function [v,t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function obj = step(obj) + obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, obj.F, obj.coeffs); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file
--- a/+time/Rungekutta4.m Wed Sep 19 16:32:05 2018 +0200 +++ b/+time/Rungekutta4.m Thu Sep 20 17:51:19 2018 +0200 @@ -39,7 +39,7 @@ end function obj = step(obj) - obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F); + 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
--- a/+time/Rungekutta4RV.m Wed Sep 19 16:32:05 2018 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,50 +0,0 @@ -classdef Rungekutta4RV < time.Timestepper - properties - F - k - t - v - m - n - - % Additional members used for the RV update - RV - end - - - methods - function obj = Rungekutta4RV(F, k, t0, v0, RV) - obj.F = F; - obj.k = k; - obj.t = t0; - obj.v = v0; - obj.m = length(v0); - obj.n = 0; - obj.RV = RV; - end - - function [v, t] = getV(obj) - v = obj.v; - t = obj.t; - end - - function state = getState(obj) - [residual, u_t, f_x] = obj.RV.getResidual(); - state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'f_x', f_x, 'viscosity', obj.RV.getViscosity(), 't', obj.t); - end - - function obj = step(obj) - obj.v = time.rk4.rungekutta_6RV(obj.v, obj.t, obj.k, obj.F, obj.RV); - 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/Rungekutta4SecondOrder.m Wed Sep 19 16:32:05 2018 +0200 +++ b/+time/Rungekutta4SecondOrder.m Thu Sep 20 17:51:19 2018 +0200 @@ -99,7 +99,7 @@ end function obj = step(obj) - obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); + obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F); obj.t = obj.t + obj.k; obj.n = obj.n + 1; end
--- a/+time/Rungekutta4proper.m Wed Sep 19 16:32:05 2018 +0200 +++ b/+time/Rungekutta4proper.m Thu Sep 20 17:51:19 2018 +0200 @@ -26,7 +26,7 @@ end function obj = step(obj) - obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F); + 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/RungekuttaRV.m Thu Sep 20 17:51:19 2018 +0200 @@ -0,0 +1,50 @@ +classdef RungekuttaRV < time.Timestepper + properties + F % RHS of the ODE + k % Time step + t % Time point + v % Solution vector + n % Time level + RV % Residual Viscosity + coeffs % The coefficents used for the RK time integration + end + + methods + function obj = RungekuttaRV(F, k, t0, v0, RV, order) + obj.F = F; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + obj.RV = RV; + % Extract the coefficients for the specified order + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = time.rk.butcherTableau(order); + obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); + end + + function [v, t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function state = getState(obj) + [residual, u_t, f_x] = obj.RV.getResidual(); + state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'f_x', f_x, 'viscosity', obj.RV.getViscosity(), 't', obj.t); + end + + function obj = step(obj) + obj.v = time.rk.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.coeffs); + 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