Mercurial > repos > public > sbplib

changeset 887:50d5a3843099 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename package rk4 to rk
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 15 Nov 2018 16:42:58 -0800
parents 8894e9c49e40
children 8732d6bd9890
files +time/+rk/get_rk4_time_step.m +time/+rk/rk4_stability.m +time/+rk/rungekutta_4.m +time/+rk/rungekutta_6.m +time/+rk4/get_rk4_time_step.m +time/+rk4/rk4_stability.m +time/+rk4/rungekutta_4.m +time/+rk4/rungekutta_6.m +time/Rk4SecondOrderNonlin.m +time/Rungekutta4.m +time/Rungekutta4SecondOrder.m +time/Rungekutta4proper.m
diffstat 12 files changed, 124 insertions(+), 124 deletions(-) [+]
line wrap: on
line diff
diff -r 8894e9c49e40 -r 50d5a3843099 +time/+rk/get_rk4_time_step.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/get_rk4_time_step.m	Thu Nov 15 16:42:58 2018 -0800
@@ -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
diff -r 8894e9c49e40 -r 50d5a3843099 +time/+rk/rk4_stability.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/rk4_stability.m	Thu Nov 15 16:42:58 2018 -0800
@@ -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
diff -r 8894e9c49e40 -r 50d5a3843099 +time/+rk/rungekutta_4.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/rungekutta_4.m	Thu Nov 15 16:42:58 2018 -0800
@@ -0,0 +1,10 @@
+% 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
diff -r 8894e9c49e40 -r 50d5a3843099 +time/+rk/rungekutta_6.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/rungekutta_6.m	Thu Nov 15 16:42:58 2018 -0800
@@ -0,0 +1,31 @@
+% 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 8894e9c49e40 -r 50d5a3843099 +time/+rk4/get_rk4_time_step.m
--- a/+time/+rk4/get_rk4_time_step.m	Thu Nov 15 16:36:21 2018 -0800
+++ /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
diff -r 8894e9c49e40 -r 50d5a3843099 +time/+rk4/rk4_stability.m
--- a/+time/+rk4/rk4_stability.m	Thu Nov 15 16:36:21 2018 -0800
+++ /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
diff -r 8894e9c49e40 -r 50d5a3843099 +time/+rk4/rungekutta_4.m
--- a/+time/+rk4/rungekutta_4.m	Thu Nov 15 16:36:21 2018 -0800
+++ /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
diff -r 8894e9c49e40 -r 50d5a3843099 +time/+rk4/rungekutta_6.m
--- a/+time/+rk4/rungekutta_6.m	Thu Nov 15 16:36:21 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 8894e9c49e40 -r 50d5a3843099 +time/Rk4SecondOrderNonlin.m
--- a/+time/Rk4SecondOrderNonlin.m	Thu Nov 15 16:36:21 2018 -0800
+++ b/+time/Rk4SecondOrderNonlin.m	Thu Nov 15 16:42:58 2018 -0800
@@ -67,7 +67,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
diff -r 8894e9c49e40 -r 50d5a3843099 +time/Rungekutta4.m
--- a/+time/Rungekutta4.m	Thu Nov 15 16:36:21 2018 -0800
+++ b/+time/Rungekutta4.m	Thu Nov 15 16:42:58 2018 -0800
@@ -43,7 +43,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
diff -r 8894e9c49e40 -r 50d5a3843099 +time/Rungekutta4SecondOrder.m
--- a/+time/Rungekutta4SecondOrder.m	Thu Nov 15 16:36:21 2018 -0800
+++ b/+time/Rungekutta4SecondOrder.m	Thu Nov 15 16:42:58 2018 -0800
@@ -105,7 +105,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
diff -r 8894e9c49e40 -r 50d5a3843099 +time/Rungekutta4proper.m
--- a/+time/Rungekutta4proper.m	Thu Nov 15 16:36:21 2018 -0800
+++ b/+time/Rungekutta4proper.m	Thu Nov 15 16:42:58 2018 -0800
@@ -31,7 +31,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