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.

请先登录，再进行评论。

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

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

采纳的回答

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

更多回答（1 个）

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

另请参阅

类别

标签

Community Treasure Hunt

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

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

采纳的回答

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

更多回答（1 个）

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

另请参阅

类别

标签

Community Treasure Hunt

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

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

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