What Are the Correct Steps to Plot Stability Regions in Numerical Methods?

[Ally1]

In many cases one finds code such as
x = linspace(-3,1.5,200);
y = linspace(-3.5,3.5,200);
[X,Y] = meshgrid(x,y);
Z = X +Y*i;
%Euler's Method
M = abs(1+Z);
[c,h] = contour(X,Y,M,[1,1]);
set(h,'linewidth',2,'edgecolor','g')
to plot the stability region of the Euler's Method, where in fact the definition of is Z = lambda*h, where h is a step-size and lambda an eigenvalue.

1. How do I modify this program the above in terms of Z = lambda*h?

2. How do I plot different regions on same plot?

Thanking you in advance.

 

[Dethrone]

Hi Ally,

For 2), you can plot different/multiple regions on the same plot by using the "hold on" and "hold off" command.

 

[Ally1]

Ok thanks but I am more concern about point 1).

 

[Ackbach]

Ok thanks but I am more concern about point 1).

Well, if $\lambda$ is an eigenvalue, then it's the eigenvalue of some operator. I have one possibility in mind for what that operator is, but what do you think it is?

 

[Ally1]

Thank for your respond.

To clear my question, let me take an example of the Exponential Time Differencing-Runge-Kutta Method [F. de la Hoz and F. Vadillo Computer Physics Communications 179 (2008) 449–456 ], I tried to code the stability regiong of this method as you can see below but the code is not giving me the correct regions within each other.
function stability_of_cancer_hiv(y)
%
ROOT = []; ROOT0 = []; lambda = 0.05; h = 0.6; z = lambda*h;
%
if y==0;
% y=-0.5*h;
c0 = exp(y);
c1 = 1 + y + 1/2*y^2 + 1/6*y^3 + 13/320*y^4 + 7/960*y^5;
c2 = 1/2 + 1/2*y + 1/4*y^2 + 247/2880*y^3 + 131/5760*y^4 + 479/96768*y^5;
c3=1/6+1/6*y+61/720*y^2+1/36*y^3+1441/241920*y^4+67/120960*y^5;
c4=1/24+1/32*y+7/640*y^2+19/11520*y^3-25/64512*y^4-311/860160*y^5;
else
c0=exp(y);
c1=-4/y^3+8*exp(y/2)/y^3-8*exp(3*y/2)/y^3+4*exp(2*y)/y^3-1/y^2+4*exp(y/2)/y^2-...
6*exp(y)/y^2+4*exp(3*y/2)/y^2-exp(2*y)/y^2;
c2=-8/y^4+16*exp(y/2)/y^4-16*exp(3*y/2)/y^4+8*exp(2*y)/y^4-5/y^3+12*exp(y/2)/y^3-...
10*exp(y)/y^3+4*exp(3*y/2)/y^3-exp(2*y)/y^3-1/y^2+4*exp(y/2)/y^2-3*exp(y)/y^2;
c3=4/y^5-16*exp(y/2)/y^5+16*exp(y)/y^5+8*exp(3*y/2)/y^5-20*exp(2*y)/y^5+8*exp(5*y/2)/y^5+...
2/y^4-10*exp(y/2)/y^4+16*exp(y)/y^4-12*exp(3*y/2)/y^4+6*exp(2*y)/y^4-2*exp(5*y/2)/y^4+...
4*exp(y)/y^3-2*exp(3*y/2)/y^3;
c4=8/y^6-24*exp(y/2)/y^6+16*exp(y)/y^6+16*exp(3*y/2)/y^6-24*exp(2*y)/y^6+8*exp(5*y/2)/y^6+...
6/y^5-18*exp(y/2)/y^5+20*exp(y)/y^5-12*exp(3*y/2)/y^5+6*exp(2*y)/y^5-2*exp(5*y/2)/y^5+...
2/y^4-6*exp(y/2)/y^4+6*exp(y)/y^4-2*exp(3*y/2)/y^4;
end
for theta = 1:200;
p0 = [z^4/24 z^3/6 z^2/2 z (1-exp(1i*theta))];
ROOT1 = roots(p0);
ROOT = [ROOT ROOT1];
%
p1 = [c4*z^4/24 c3*z^4/6 c2*z^4/2 c1*z c0*(1-exp(1i*theta))];
ROOT2 = roots(p1);
ROOT0 = [ROOT0 ROOT2];
end
% disp(ROOT);
figure(1); plot(real(ROOT),imag(ROOT),'k','linewidth',0.1);
figure(2); plot(real(ROOT0),imag(ROOT0),'k','linewidth',1);