Stability boundaries - intersection with real and Imaginary axis

I am plotting the stability regions (please see the code below).

How can I calculate the intersections with real and imaginary axis for each stability region, please?

I was trying to filter the data in the matrix M (size 200x200), but unfortunately, I was not very successulf:

For example, I was trying to filter the matrix data to get the intersection with real axis:

I tried to get the entries where imaginary part of M equals to zero, and then get the maximum real value, but it seems it does not work correclty:

zeroImag = M(imag(abs(M)) == 0);
maxReal = max(real(abs(zeroImag)))

The code for plotting the stability regions:

x = linspace(-5, 2, 200);
y = linspace(-3.5, 3.5, 200);
[X,Y] = meshgrid(x,y);
Z = X + 1i*Y;

% Stability function 
% R(z) = 1 + z + z^2/2! + z^3/3! + z^4/4 + ...
nom = Z;
denom = 1;
M = 1;

figure;
color = {'b','r','g','m','c','k'};

for j=2:5
    M = M + nom/denom;
    nom = Z.^j; 
    denom = j * denom;
    % plot the stability region of order "j" 
    contour(x,y,abs(M),[1 1], 'Color', color{j}, 'LineWidth',2);
    hold on;  
end
grid on;

% plot real and imaginary axis
plot([-5, 2],[0 0],'k-')
hold on;
plot([0 0],[-3.5, 3.5],'k-')

xlabel('Re')
ylabel('Im')
title('Runge-Kutta - orders 1,2,3,4')

Answer 1

[Torsten]

You can determine the roots analytically because the polynomials have order <= 4.

Here is a numerical solution for imag(Z) = 0. The solution for real(Z)=0 can be obtained similarly.

format long
M2 = @(x,y) 1 + (x + 1i*y);
M3 = @(x,y) 1 + (x + 1i*y) + (x + 1i*y)^2/2;
M4 = @(x,y) 1 + (x + 1i*y) + (x + 1i*y)^2/2 + (x + 1i*y)^3/6;
M5 = @(x,y) 1 + (x + 1i*y) + (x + 1i*y)^2/2 + (x + 1i*y)^3/6 + (x + 1i*y)^4/24;
options = optimset('TolX',1e-12,'TolF',1e-12,'Display','none');
fun = @(x) abs(M2(x,0))-1;
sol = fsolve(fun,-3,options)

sol = 
    -2

fun(sol)

ans = 
     0

fun = @(x) abs(M3(x,0))-1;
sol = fsolve(fun,-3,options)

sol = 
  -2.000000000000003

fun(sol)

ans = 
     2.664535259100376e-15

fun = @(x) abs(M4(x,0))-1;
sol = fsolve(fun,-3,options)

sol = 
  -2.512745326618329

fun(sol)

ans = 
    -4.440892098500626e-16

fun = @(x) abs(M5(x,0))-1;
sol = fsolve(fun,-3,options)

sol = 
  -2.785293563405308

fun(sol)

ans = 
     4.041211809635570e-14

[Carci]

Torsten, thank you very much for your answer!

And if I need to evaluate higher-orders (for example the 10th order), then it is probably better to check the values in the plot manually using the"Brush/Select data" tool , is that correct?

[Torsten]

I'd use "roots" for the real or imaginary part of the polynomial abs(M).^2 - 1 to get the intersections with the real or imaginary axis.

Example:

syms x y
assume (x,'real')
assume (y,'real')
order = 10;
M = sum((x+1i*y).^(0:order)./factorial(0:order));
p = M*M' - 1;
p_real = sym2poly(subs(p,y,0));
r_real = roots(p_real);
r_real = unique(r_real(abs(imag(r_real))<eps))

r_real = 2×1

-5.0695 0

p_imag = sym2poly(subs(p,x,0));
r_imag = roots(p_imag);
r_imag = unique(r_imag(abs(imag(r_imag))<eps))

r_imag = 5×1

-5.2619 -3.4324 0 3.4324 5.2619

[Carci]

Torsten, thank you very much again! :-) It really helped me a lot.