Mercurial > repos > public > sbplib
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
--- 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
--- 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