Stability boundaries - intersection with real and Imaginary axis

12 次查看（过去 30 天）
Carci
Carci 2022-9-10
评论： Carci 2022-9-11
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')

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-9-10
编辑：Torsten 2022-9-10
Ran in:
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
  3 个评论
显示 1更早的评论
Torsten
Torsten 2022-9-10
编辑：Torsten 2022-9-10
Ran in:
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by