Plotting stability regions for Euler's method with Z = lambda*h

  • Context:
  • Thread starter Thread starter Ally1
  • Start date Start date
  • Tags Tags
    Plotting Stability
Join the discussion
Registration is free. Ask a follow-up in this thread, or start your own.
Sign up Log in
4 replies · 4K views
Ally1
Messages
3
Reaction score
0
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.
 

Discussion

Physics news on Phys.org
Hi Ally,

For 2), you can plot different/multiple regions on the same plot by using the "hold on" and "hold off" command.
 
Ok thanks but I am more concern about point 1).
 
Ally said:
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?
 
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);