Figure 2.28:

Stability regions for Runge-Kutta methods; \dot {x}=\lambda x.

Code for Figure 2.28

Text of the GNU GPL.

main.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
% Program to plot stability plots for Runge Kutta methods
% For details on plotting stability regions for One-step Method,
% see pg 163 of Randall J. LeVeque book on " Finite Difference
% Methods for Ordinary and Partial Differential Equations


% In case on one-step method, any time-step n+1 can be written as
% u(n+1) = R(z)*u(n)
% where R(z) is some function of z = lambda*dt

% The region of absolute stability for a one-step method is simply the
% region where |R(z)| <= 1


% We demonstrate here stability curve for RK2, RK3 and RK4 methods.

% Setting up the grid points

axisbox = [-4 4 -4 4];
xa = axisbox(1); xb = axisbox(2); ya = axisbox(3); yb = axisbox(4);
nptsx = 51;
nptsy = 51;
x = linspace(xa,xb,nptsx);
y = linspace(ya,yb,nptsy);
[X,Y] = meshgrid(x,y);
Z = X + i*Y;

% 2nd order Runge-Kutta method
RK2 = @(z) (1 + z + 1/2*z.^2);

% 3rd order Runge-Kutta method
RK3 = @(z) (1 + z + 1/2*z.^2 + 1/6*z.^3);

% 4th order Runge-Kutta method
RK4 = @(z) (1 + z + 1/2*z.^2 + 1/6*z.^3 + 1/24*z.^4);

% Evaluation of R(z) inside axisbox
Rval2 = RK2(Z); Rval3 = RK3(Z); Rval4 = RK4(Z);
%----------------------------------------------------------------------
% evaluation of |R(z)| for finding S, region of absolute stability:

Rabs2 = abs(Rval2); Rabs3 = abs(Rval3); Rabs4 = abs(Rval4);

% Plotting the stability curve using the contour plot of MATLAB

figure()
clf()
hold on
% plot axes:
plot([xa xb],[0 0],'k-','LineWidth',2)
plot([0 0],[ya yb],'k-','LineWidth',2)
contour(x,y,Rabs2,[1 1],'r-','LineWidth',2);
contour(x,y,Rabs3,[1 1],'b-','LineWidth',2);
contour(x,y,Rabs4,[1 1],'g-','LineWidth',2);
legend('RK2','RK3','RK4')
title('Region of absolute stability','FontSize',15)
axis([xa xb ya yb])
set(gca,'FontSize',15)
hold off

% store data
c2 = contourc(x,y,Rabs2,[1 1])';
c2(1,:) = [];
c3 = contourc(x,y,Rabs3,[1 1])';
c3(1,:) = [];
c4 = contourc(x,y,Rabs4,[1 1])';
c4(1,:) = [];
save "RK_stability.dat" c2 c3 c4