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

62 ビュー (過去 30 日間)
KS 2024 年 5 月 18 日
コメント済み: Paul 2024 年 5 月 21 日 12:32
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.

サインインしてコメントする。

### 採用された回答

Walter Roberson 2024 年 5 月 19 日

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
##### 7 件のコメント5 件の古いコメントを表示5 件の古いコメントを非表示
KS 2024 年 5 月 21 日 11:25
@Walter Roberson Many thanks
Paul 2024 年 5 月 21 日 12:32
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).

サインインしてコメントする。

### その他の回答 (1 件)

Image Analyst 2024 年 5 月 19 日
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 件の古いコメントを非表示
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 日 16:02
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.

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeData Distribution Plots についてさらに検索

### Community Treasure Hunt

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

Start Hunting!

Translated by