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

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

[John D'Errico]

All you show is a picture. No code. What instability? How is this even a question about MATLAB?

If you insist, I'll try the MATLAB crystal ball. Yes, it tells me something now. I think it says you will inherit some money soon. Well, maybe your heirs will inherit some money. Darn foggy crystal balls. It should teach me a lesson about buying a discount imported crystal ball online.

[Image Analyst]

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]

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

Answer 1

[Walter Roberson]

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

[Paul]

I tried a nearly exact same approach and it didn't work for me. On quick glance, the only difference in you answer and my attempt are how you treated A0 and A1. I had

A0  = Q(1.5207);
A1  = Q(0.853721-1.816893i);

Any idea why my A0 and A1 didn't work? Or why did you implement A0 and A1 as you did?

[KS]

@Walter Roberson thank you for your support. Though it solved the problem of unwanted data jumps, I still haven't got the expected result. What step can be taken next?

[KS]

Hi @Paul , I tried your approach and it works. However, you may try with quotes. For example,

A0 = Q('1.5207');

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

The symbolic (sym) numbers are represented in exact rational form (gives exact calculations) rather than decimal approxmiations that in the case of floating-point numbers.

[Sam Chak]

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

[KS]

@Sam Chak Many thanks for your support. It resolved my problem

[KS]

@Walter Roberson Many thanks

[Paul]

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]

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]

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-');

[KS]

@Image Analyst thanks, I know that. My question is like - the expected Plot is theoretically supposed to be like "green color," but my code ("pink color") has got some unwanted jumps up when theta>40, which is not supposed to be.

[Sam Chak]

@KS, The problem seems to have started here.

A0  = 1.5207;
A1  = 0.853721-1.816893i;
A2  = 1;
th  = 40:0.001:42;
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)); 
plot(th, th2), grid, xlabel \theta, ylabel th2

Warning: Imaginary parts of complex X and/or Y arguments ignored.

c2  = cos(th2);
plot(th, c2), grid, xlabel \theta, ylabel c2

Warning: Imaginary parts of complex X and/or Y arguments ignored.

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

[KS]

@Sam Chak thanks for details. Exactly, I've also noticed. How to solve this problem?

[Sam Chak]

@KS, This term "(A1/A2).*sin(th1)" is the culprit.

A0  = 1.5207;
A1  = 0.853721-1.816893i;
A2  = 1;
th  = 40:0.001:42;
th0 = th*pi/180;    % deg to rad
c0  = cos(th0);
th1 = asin((A0/A1).*sin(th0));
c1  = cos(th1);
plot(th, (A1/A2).*sin(th1)), grid, xlabel \theta

Warning: Imaginary parts of complex X and/or Y arguments ignored.

Let me quote @John D'Errico: "There is NO REAL angle that returns a sin of a number such that the sin is greater than 1, or less than minus 1. As such, asin can NEVER return a real result for inputs outside of that range."

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

[Sam Chak]

@KS, The green curve after 41.1165° seems to be the lower envelope of P. Can you determine the values of this lower envelope?

[Image Analyst]

@KS I'm not sure what the "problem" is. It's plotting what you've asked it to plot (in blue or green, if asked, not pink like you said). What are you expecting or what would you like the plot to look like if there was no "problem"? Maybe you should just plot the real or imaginary component alone.

[KS]

@Sam Chak How to determine the values of lower enevelope?

[Image Analyst]

Just curious, what do those values actually mean, like in a real world physical sense? And why do you need them?

[Sam Chak]

@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]

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

[Sam Chak]

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.