フィルターのクリア

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

62 ビュー (過去 30 日間)
KS
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 件のコメント
Image Analyst
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
KS 2024 年 5 月 19 日
Thank you @John D'Errico and @Image Analyst. I've updated the question with code.

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

採用された回答

Walter Roberson
Walter Roberson 2024 年 5 月 19 日
移動済み: 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 件のコメント
KS
KS 2024 年 5 月 21 日 11:25
@Walter Roberson Many thanks
Paul
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
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 件のコメント
KS
KS 2024 年 5 月 20 日
@Sam Chak thanks for your suggesstion. I've tried with "envelope" function, doesn't work in this case.
Sam Chak
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