changeset 994:2f89959fb9f0 feature/timesteppers

Implement method to get gain from butcher tableu
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 09 Jan 2019 12:14:30 +0100
parents 44e7e497c3b7
children 10c5eda235b7
files +time/+rk/ButcherTableau.m +time/+rk/rk4_stability.m
diffstat 2 files changed, 12 insertions(+), 58 deletions(-) [+]
line wrap: on
line diff
diff -r 44e7e497c3b7 -r 2f89959fb9f0 +time/+rk/ButcherTableau.m
--- a/+time/+rk/ButcherTableau.m	Wed Jan 09 11:14:16 2019 +0100
+++ b/+time/+rk/ButcherTableau.m	Wed Jan 09 12:14:30 2019 +0100
@@ -27,6 +27,18 @@
             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
 
diff -r 44e7e497c3b7 -r 2f89959fb9f0 +time/+rk/rk4_stability.m
--- a/+time/+rk/rk4_stability.m	Wed Jan 09 11:14:16 2019 +0100
+++ /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