Mercurial > repos > public > sbplib

comparison +time/+rk4/rk4_stability.m @ 0:48b6fb693025

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial commit.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 17 Sep 2015 10:12:50 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:48b6fb693025
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