7.5. Example MATLAB Programs#
7.5.1. Script 1#
%6G6Z3002 - Computational methods in ODEs
%Interval of Absolute Stability (script 1)
%Matlab script for producing a table for interval of absolute stability
%using Root Locus Method of a given LMS method- here AM2-step
%Note: the result is printed in columns, and not in rows as given in the
%lecture notes
out=[];
for h=1:-1:-7
p=[(1-5*h/12) -(1+8*h/12) (h/12)];
r=abs(roots(p));
out=[h;r(1);r(2)];
fprintf('%6.1f %6.2f %6.2f\n',out)
end
Output:
1.0 2.81 0.05
0.0 0.00 1.00
-1.0 0.39 0.15
-2.0 0.41 0.22
-3.0 0.62 0.18
-4.0 0.78 0.16
-5.0 0.91 0.15
-6.0 1.00 0.14
-7.0 1.07 0.14
7.5.2. Script 2#
%6G6Z3002 - Computational methods in ODEs
%Interval of Absolute Stability (script 2)
%Matlab script for producing a table for interval of absolute stability
%using Root Locus Method of a given LMS method- here AM2-step
%Note: the result is printed in rows as given in the lecture notes
%
out=[];
Hlam=[];
r1=[];
r2=[];
count=0;
h=2;
while h >= -7
count=count+1;
h=h-1;
p=[(1-5*h/12) -(1+8*h/12) (h/12)];
r=roots(p);
Hlam(count)=h;
r1(count)=r(1);
r2(count)=r(2);
end
disp(count-1);
for i=count-1:-1:1
fprintf('%5.2f\t',Hlam(i))
end
fprintf('\n')
fprintf('----------------------------------------------------------\n');
for i=count-1:-1:1
fprintf('%5.2f\t',abs(r1(i)))
end
fprintf('\n')
for i=count-1:-1:1
fprintf('%5.2f\t',abs(r2(i)))
end
fprintf('\n')
Output
-7.00 -6.00 -5.00 -4.00 -3.00 -2.00 -1.00 0.00 1.00
-------------------------------------------------------------
1.07 1.00 0.91 0.78 0.62 0.41 0.39 0.00 2.81
0.14 0.14 0.15 0.16 0.18 0.22 0.15 1.00 0.05
7.5.3. Script 3#
%6G6Z3002 - Computational methods in ODEs
%Plot of Region of Absolute Stability (script 1) AB 2step
clc
clear
syms theta
hlamda = (2*(exp(2i*theta) - exp(1i*theta)))/(3* exp(1i*theta)-1);
H=[];
THETA=[];
for i = 0:pi/30:2*pi
THETA = [THETA i];
H = [H subs(hlamda,theta,i)];
end
a = max(size(THETA));
fprintf(' theta real(hlamda) imag(hlamda)\n')
fprintf('-----------------------------------------------\n')
for i = 1:a
fprintf('%6i %8.2f %12.2f\n',round((THETA(i)/pi)*180),double(real(H(i))),double(imag(H(i))))
end
figure(1)
plot(real(H),imag(H),'LineWidth',2);
set(gca,'FontSize',12);
xlabel('Re (h\lambda)');
ylabel('Im (h\lambda)');
title('Region of Absolute Stability for 2-step AB method');
axis([-1 0.4 -1 1]);
grid on;
Output
Fig. 7.4 Stability region of AB 2-step method.#
7.5.4. Script 4#
%A very simple version
%Region of absolute stability for the 2-step AB method
w=exp(1i*linspace(0,2*pi));
z=2*(w.^2-w)./(3*w-1);
plot(z)
7.5.5. Script 5#
%Plot of Region of Absolute Stability (script 2)-AB_2step
clc
clear
syms theta
z=exp(1i*theta);
hlamda = (2*(z^2 - z))/(3*z-1);
H=[];
THETA=[];
for i = 0:pi/30:2*pi
THETA = [THETA i];
H = [H subs(hlamda,theta,i)];
end
a = max(size(THETA));
fprintf(' \x03B8 Re(h\x03bb) Im(h\x03bb)\n')
fprintf('------------------------------\n')
for i = 1:a
fprintf('%4i %8.2f %10.2f\n', ...
round((THETA(i)/pi)*180),double(real(H(i))),double(imag(H(i))))
end
fill(real(H),imag(H),'c');
set(gca,'FontSize',12);
xlabel('Re (h\lambda)');
ylabel('Im (h\lambda)');
title('Region of Absolute Stability for 2-step AB method');
axis([-1 0.4 -1 1]);
grid on;
Output
θ Re(hλ) Im(hλ)
------------------------------
0 0.00 0.00
6 -0.00 0.10
12 -0.00 0.21
18 -0.00 0.30
24 -0.01 0.39
30 -0.01 0.47
36 -0.03 0.54
42 -0.05 0.61
. . .
. . .
. . .
168 -0.99 0.16
174 -1.00 0.08
180 -1.00 0.00
. . .
. . .
. . .
276 -0.34 -0.80
282 -0.29 -0.80
288 -0.23 -0.79
294 -0.19 -0.77
300 -0.14 -0.74
306 -0.10 -0.71
312 -0.07 -0.66
318 -0.05 -0.61
324 -0.03 -0.54
330 -0.01 -0.47
336 -0.01 -0.39
342 -0.00 -0.30
348 -0.00 -0.21
354 -0.00 -0.10
360 0.00 0.00
Fig. 7.5 Stability region of AB 2-step method.#