|
0
|
1 function rk_stability()
|
|
|
2 ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4));
|
|
|
3 circ = @(z)(abs(z));
|
|
|
4
|
|
|
5
|
|
|
6 % contour(X,Y,z)
|
|
|
7 ax = [-4 2 -3 3];
|
|
|
8 % hold on
|
|
|
9 fcontour(ruku4,[1,1],[-3, 0.6],[-3.2, 3.2])
|
|
|
10 hold on
|
|
|
11 r = 2.6;
|
|
|
12 fcontour(circ,[r,r],[-3, 0.6],[-3.2, 3.2],'r')
|
|
|
13 hold off
|
|
|
14 % contour(X,Y,z,[1,1],'b')
|
|
|
15 axis(ax)
|
|
|
16 title('4th order Runge-Kutta stability region')
|
|
|
17 xlabel('Re')
|
|
|
18 ylabel('Im')
|
|
|
19 axis equal
|
|
|
20 grid on
|
|
|
21 box on
|
|
|
22 hold off
|
|
|
23 % surf(X,Y,z)
|
|
|
24
|
|
|
25
|
|
|
26 rk4roots()
|
|
|
27 end
|
|
|
28
|
|
|
29 function fcontour(f,levels,x_lim,y_lim,opt)
|
|
|
30 default_arg('opt','b')
|
|
|
31 x = linspace(x_lim(1),x_lim(2));
|
|
|
32 y = linspace(y_lim(1),y_lim(2));
|
|
|
33 [X,Y] = meshgrid(x,y);
|
|
|
34 mu = X+ 1i*Y;
|
|
|
35
|
|
|
36 z = f(mu);
|
|
|
37
|
|
|
38 contour(X,Y,z,levels,opt)
|
|
|
39
|
|
|
40 end
|
|
|
41
|
|
|
42
|
|
|
43 function rk4roots()
|
|
|
44 ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4));
|
|
|
45 % Roots for real evalues:
|
|
|
46 F = @(x)(abs(ruku4(x))-1);
|
|
|
47 real_x = fzero(F,-3);
|
|
|
48
|
|
|
49 % Roots for imaginary evalues:
|
|
|
50 F = @(x)(abs(ruku4(1i*x))-1);
|
|
|
51 imag_x1 = fzero(F,-3);
|
|
|
52 imag_x2 = fzero(F,3);
|
|
|
53
|
|
|
54
|
|
|
55 fprintf('Real x = %f\n',real_x)
|
|
|
56 fprintf('Imag x = %f\n',imag_x1)
|
|
|
57 fprintf('Imag x = %f\n',imag_x2)
|
|
|
58 end |
