What are the possible reasons for data jumps in the plot?

Question

KS 2024-5-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2120266-what-are-the-possible-reasons-for-data-jumps-in-the-plot

编辑： KS 2024-6-24

Plot_Data Jump.png

What might be the cause of such numerical instability and how can it be solved? Note: Expected plot should be as like "Green" color.

A0 = 1.5207;

A1 = 0.853721-1.816893i;

A2 = 1;

th = 0:0.5:90;

th0 = th*pi/180; % deg to rad

c0 = cos(th0);

th1 = asin((A0/A1).*sin(th0));

c1 = cos(th1);

th2 = asin((A1/A2).*sin(th1));

c2 = cos(th2);

b = 2.*pi.*50.*A1.*c1./500;

M1 = (A1.*c0 - A0.*c1)./(A1.*c0 + A0.*c1);

M2 = (A2.*c1 - A1.*c2)./(A2.*c1 + A1.*c2);

M = (M1 + (M2.*exp(-2i.*b)))./(1+(M1.*M2.*exp(-2i.*b)));

P = abs(M).^2;

plot(th, P), grid, xlabel \theta, ylabel P

3 个评论
显示 1更早的评论隐藏 1更早的评论

Image Analyst 2024-5-18

Yeah @John D'Errico I think the Crystal Ball Toolbox is the wrong one here. I think you need the Mind Reading Toolbox. Unfortunately it's not released yet and only @Walter Roberson has an early alpha release of it.

KS 2024-5-19

Thank you @John D'Errico and @Image Analyst. I've updated the question with code.

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2024-5-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2120266-what-are-the-possible-reasons-for-data-jumps-in-the-plot#answer_1460001

移动：Walter Roberson 2024-5-19

在 MATLAB Online 中打开

Q = @(v) sym(v);

Pi = Q(pi);

A0 = Q(15207)/Q(10)^4;

A1 = Q(0853721)/Q(10)^6-Q(1816893)/Q(10)^6*1i;

A2 = Q(1);

th = Q(0:0.5:90);

th0 = th*Pi/180; % deg to rad

c0 = cos(th0);

th1 = asin((A0/A1).*sin(th0));

c1 = cos(th1);

th2 = asin((A1/A2).*sin(th1));

c2 = cos(th2);

b = 2.*Pi.*50.*A1.*c1./500;

M1 = (A1.*c0 - A0.*c1)./(A1.*c0 + A0.*c1);

M2 = (A2.*c1 - A1.*c2)./(A2.*c1 + A1.*c2);

M = (M1 + (M2.*exp(-2i.*b)))./(1+(M1.*M2.*exp(-2i.*b)));

P = abs(M).^2;

plot(th, P), grid, xlabel \theta, ylabel P

8 个评论
显示 6更早的评论隐藏 6更早的评论

Sam Chak 2024-5-21

编辑：Sam Chak 2024-5-21

在 MATLAB Online 中打开

Great! Using @KS–@Paul's solution, it produces the green curve (lower envelope) as desired.

A0 = 1.5207;

A1 = 0.853721-1.816893i;

A2 = 1;

th = 0:0.05:90;

th0 = th*pi/180; % deg to rad

c0 = cos(th0);

th1 = asin((A0/A1).*sin(th0));

c1 = cos(th1);

th2 = asin((A1/A2).*sin(th1));

c2 = cos(th2);

b = 2.*pi.*50.*A1.*c1./500;

M1 = (A1.*c0 - A0.*c1)./(A1.*c0 + A0.*c1);

M2 = (A2.*c1 - A1.*c2)./(A2.*c1 + A1.*c2);

M = (M1 + (M2.*exp(-2i.*b)))./(1+(M1.*M2.*exp(-2i.*b)));

P = abs(M).^2;

plot(th, P, 'color', '#F09494'), grid, xlabel \theta, ylabel P

hold on

Q = @(v) sym(v);

Pi = Q(pi);

% A0 = Q(15207)/Q(10)^4;

% A1 = Q(0853721)/Q(10)^6 - Q(1816893)/Q(10)^6*1i; % 0.853721 - 1.816893i;

A0 = Q('1.5207');

A1 = Q('0.853721 - 1.816893i');

A2 = Q(1);

th = Q(0:0.5:90);

th0 = th*Pi/180; % deg to rad

c0 = cos(th0);

th1 = asin((A0/A1).*sin(th0));

c1 = cos(th1);

th2 = asin((A1/A2).*sin(th1));

c2 = cos(th2);

b = 2.*Pi.*50.*A1.*c1./500;

M1 = (A1.*c0 - A0.*c1)./(A1.*c0 + A0.*c1);

M2 = (A2.*c1 - A1.*c2)./(A2.*c1 + A1.*c2);

M = (M1 + (M2.*exp(-2i.*b)))./(1+(M1.*M2.*exp(-2i.*b)));

P = abs(M).^2;

plot(th, P, 'color', '#AAE2A8', 'linewidth', 3), grid, xlabel \theta, ylabel P, grid

Paul 2024-5-21

在 MATLAB Online 中打开

Well, I guess it really wan't my solution insofar as I didn't get it to work. Maybe @Walter Roberson did something else that I missed

I see that using those single quotes defining A0 and A1 can make a small difference

Q = @(v) sym(v);

Q('1.5207') - Q(1.5207)

ans =

0.0

Q('0.853721 - 1.816893i') - Q(0.853721 - 1.816893i)

ans =

Is that the reason this solution covers the lower part of the envelope and Walter's covers the upper part? Seems very sensitive.

Another approach is to form P entirely symbolically, and the substitute the nuerical constants at the end, which sometimes takes advantage of opportunities to simplify expressions that might not otherwise be seen with a bunch of VPA numbers floating around.

Pi = Q(pi);

syms A0 A1

% A0 = Q(15207)/Q(10)^4;

% A1 = Q(0853721)/Q(10)^6 - Q(1816893)/Q(10)^6*1i; % 0.853721 - 1.816893i;

%A0 = Q('1.5207');

%A1 = Q('0.853721 - 1.816893i');

A2 = Q(1);

th = Q(0:0.5:90);

%th0 = th*Pi/180; % deg to rad

syms th0

c0 = cos(th0);

th1 = asin((A0/A1).*sin(th0));

c1 = cos(th1);

th2 = asin((A1/A2).*sin(th1));

c2 = cos(th2);

b = 2.*Pi.*50.*A1.*c1./Q(500);

M1 = (A1.*c0 - A0.*c1)./(A1.*c0 + A0.*c1);

M2 = (A2.*c1 - A1.*c2)./(A2.*c1 + A1.*c2);

M = (M1 + (M2.*exp(-2i.*b)))./(1+(M1.*M2.*exp(-2i.*b)));

M = simplify(M);

P = abs(M).^2;

P = subs(P,[A0 A1],[Q('1.5207'), Q('0.853721 - 1.816893i')]);

symvar(P)

ans =

P = subs(P,th0,th*Pi/180);

plot(th, P, 'color', '#AAE2A8', 'linewidth', 3), grid, xlabel \theta, ylabel P, grid

And now we have a different result (I don't know why the grid command doesn't work).

KS 2024-6-24

编辑：KS 2024-6-24

ref.png

How to read above "A0" and "A1" from spreadsheet while having more than 1 value of "A0" and "A1" [considering @Paul's explanation, attached screeenshot for quick reference]?

Note: using Quote (' ') works well in this case

请先登录，再进行评论。

Answer 2

Image Analyst 2024-5-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2120266-what-are-the-possible-reasons-for-data-jumps-in-the-plot#answer_1459951

在 MATLAB Online 中打开

I don't know where those equations come from. Is it some kind of chaos theory? Anyway, to plot in green like you asked, you need to use 'g-'

plot(th, P, 'g-');

11 个评论
显示 9更早的评论隐藏 9更早的评论

Sam Chak 2024-5-20

在 MATLAB Online 中打开

@KS, you can study @Walter Roberson's solution, which he has successfully plotted the upper envelope. Your focus and challenge now is to modify the code and plot the lower envelope.

A0 = 1.5207;

A1 = 0.853721-1.816893i;

A2 = 1;

th = 0:0.05:90;

th0 = th*pi/180; % deg to rad

c0 = cos(th0);

th1 = asin((A0/A1).*sin(th0));

c1 = cos(th1);

th2 = asin((A1/A2).*sin(th1));

c2 = cos(th2);

b = 2.*pi.*50.*A1.*c1./500;

M1 = (A1.*c0 - A0.*c1)./(A1.*c0 + A0.*c1);

M2 = (A2.*c1 - A1.*c2)./(A2.*c1 + A1.*c2);

M = (M1 + (M2.*exp(-2i.*b)))./(1+(M1.*M2.*exp(-2i.*b)));

P = abs(M).^2;

plot(th, P), grid, xlabel \theta, ylabel P

hold on

Q = @(v) sym(v);

Pi = Q(pi);

A0 = Q(15207)/Q(10)^4;

A1 = Q(0853721)/Q(10)^6 - Q(1816893)/Q(10)^6*1i; % 0.853721 - 1.816893i;

A2 = Q(1);

th = Q(0:0.5:90);

th0 = th*Pi/180; % deg to rad

c0 = cos(th0);

th1 = asin((A0/A1).*sin(th0));

c1 = cos(th1);

th2 = asin((A1/A2).*sin(th1));

c2 = cos(th2);

b = 2.*Pi.*50.*A1.*c1./500;

M1 = (A1.*c0 - A0.*c1)./(A1.*c0 + A0.*c1);

M2 = (A2.*c1 - A1.*c2)./(A2.*c1 + A1.*c2);

M = (M1 + (M2.*exp(-2i.*b)))./(1+(M1.*M2.*exp(-2i.*b)));

P = abs(M).^2;

plot(th, P, 'linewidth', 3), grid, xlabel \theta, ylabel P, grid

KS 2024-5-20

@Sam Chak thanks for your suggesstion. I've tried with "envelope" function, doesn't work in this case.

Sam Chak 2024-5-20

Hi @KS,

Could you please provide some background information on the complex-valued function you are working with? Understanding the context and purpose of the function would be helpful for us to provide accurate assistance.

请先登录，再进行评论。