polyxpoly to find intersection

i have two curves as given in the code 'p' and 'q' for diffeent values of 'x' the curves intersect at different points. Can anyone tell me why is my code not running?

clc
clf
b=linspace(0,1);
x=8;
p= x.*(sqrt(1-b)).*(((besselj(1,x.*sqrt(1-b)))./(besselj(0,x.*sqrt(1-b)))));
q = x.*(sqrt(b)).*((besselk(1,x.*sqrt(b)))./(besselk(0,x.*sqrt(b))));
plot(b,p)
hold on
plot(b,q)
[xx,yy] = polyxpoly(b,p,b,q)

Answer 1

[KSSV]

You may consider this https://in.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections?focused=5165138&tab=function to calculate the intersection points.

Save the above function in your working folder and use the below code.

clc
clf
b=linspace(0,1);
x=8;
p= x.*(sqrt(1-b)).*(((besselj(1,x.*sqrt(1-b)))./(besselj(0,x.*sqrt(1-b)))));
q = x.*(sqrt(b)).*((besselk(1,x.*sqrt(b)))./(besselk(0,x.*sqrt(b))));

L1 = [b ;p] ;
L2 = [b ;q] ; 
P = InterX(L1,L2) ; % get intersection points using the function 

plot(b,p)
hold on
plot(b,q)
plot(P(1,:),P(2,:),'*r')

[vipul kumar]

it did not work, i need the points of intersection. The numerical values of those points are to be used further.

[KSSV]

Why it will not work? Show is how did you try.

[vipul kumar]

i have given the code, if possible can you try once?

[KSSV]

I am asking the code which you tried with InterX.

[vipul kumar]

clc

clf

b=linspace(0,1);

x=2;

p= x.*(sqrt(1-b)).*(((besselj(1,x.*sqrt(1-b)))./(besselj(0,x.*sqrt(1-b)))));

q = x.*(sqrt(b)).*((besselk(1,x.*sqrt(b)))./(besselk(0,x.*sqrt(b))));

plot(b,p,b,q,P(1,:),P(2,:),'ro')

function P = InterX(L1,varargin)

%INTERX Intersection of curves

% P = INTERX(L1,L2) returns the intersection points of two curves L1

% and L2. The curves L1,L2 can be either closed or open and are described

% by two-row-matrices, where each row contains its x- and y- coordinates.

% The intersection of groups of curves (e.g. contour lines, multiply

% connected regions etc) can also be computed by separating them with a

% column of NaNs as for example

%

% L = [x11 x12 x13 ... NaN x21 x22 x23 ...;

% y11 y12 y13 ... NaN y21 y22 y23 ...]

%

% P has the same structure as L1 and L2, and its rows correspond to the

% x- and y- coordinates of the intersection points of L1 and L2. If no

% intersections are found, the returned P is empty.

%

% P = INTERX(L1) returns the self-intersection points of L1. To keep

% the code simple, the points at which the curve is tangent to itself are

% not included. P = INTERX(L1,L1) returns all the points of the curve

% together with any self-intersection points.

%

% Example:

% t = linspace(0,2*pi);

% r1 = sin(4*t)+2; x1 = r1.*cos(t); y1 = r1.*sin(t);

% r2 = sin(8*t)+2; x2 = r2.*cos(t); y2 = r2.*sin(t);

% P = InterX([x1;y1],[x2;y2]);

% plot(x1,y1,x2,y2,P(1,:),P(2,:),'ro')

% Author : NS

% Version: 3.0, 21 Sept. 2010

% Two words about the algorithm: Most of the code is self-explanatory.

% The only trick lies in the calculation of C1 and C2. To be brief, this

% is essentially the two-dimensional analog of the condition that needs

% to be satisfied by a function F(x) that has a zero in the interval

% [a,b], namely

% F(a)*F(b) <= 0

% C1 and C2 exactly do this for each segment of curves 1 and 2

% respectively. If this condition is satisfied simultaneously for two

% segments then we know that they will cross at some point.

% Each factor of the 'C' arrays is essentially a matrix containing

% the numerators of the signed distances between points of one curve

% and line segments of the other.

%...Argument checks and assignment of L2

error(nargchk(1,2,nargin));

if nargin == 1,

L2 = L1; hF = @lt; %...Avoid the inclusion of common points

else

L2 = varargin{1}; hF = @le;

end

%...Preliminary stuff

x1 = L1(1,:)'; x2 = L2(1,:);

y1 = L1(2,:)'; y2 = L2(2,:);

dx1 = diff(x1); dy1 = diff(y1);

dx2 = diff(x2); dy2 = diff(y2);

%...Determine 'signed distances'

S1 = dx1.*y1(1:end-1) - dy1.*x1(1:end-1);

S2 = dx2.*y2(1:end-1) - dy2.*x2(1:end-1);

C1 = feval(hF,D(bsxfun(@times,dx1,y2)-bsxfun(@times,dy1,x2),S1),0);

C2 = feval(hF,D((bsxfun(@times,y1,dx2)-bsxfun(@times,x1,dy2))',S2'),0)';

%...Obtain the segments where an intersection is expected

[i,j] = find(C1 & C2);

if isempty(i),P = zeros(2,0);return; end;

%...Transpose and prepare for output

i=i'; dx2=dx2'; dy2=dy2'; S2 = S2';

L = dy2(j).*dx1(i) - dy1(i).*dx2(j);

i = i(L~=0); j=j(L~=0); L=L(L~=0); %...Avoid divisions by 0

%...Solve system of eqs to get the common points

P = unique([dx2(j).*S1(i) - dx1(i).*S2(j), ...

dy2(j).*S1(i) - dy1(i).*S2(j)]./[L L],'rows')';

function u = D(x,y)

u = bsxfun(@minus,x(:,1:end-1),y).*bsxfun(@minus,x(:,2:end),y);

end

end

[vipul kumar]

here it is, please tell me if i am doing it correctly?

[KSSV]

Edited the answer.

[vipul kumar]

there is an error in line 7.

plot command line

it is showing that the variable 'P' is not defined

[KSSV]

It is running and I am able to get P. I got plot as:

P as:

0.13212 0.51661 0.63025 0.91891 0.92905

3.37471 6.23127 6.83386 8.15427 8.19651

[vipul kumar]

its not working in mine and the plot is also not correct.

Can you once send the running code here so i may check if it generates correct plot?

[KSSV]

I have already given the working code.

[vipul kumar]

ohh, you edited the parent question. My bad. Got it bruh, thanks a lot.

[vipul kumar]

There is one more thing, if you lookt at it. Since the function goes to infinity so the intersection point generated in the graph are more i.e. if you look from the right of the graph, only alternate points are taken into account. First relevent point is the rightmost and then alternate are to be taken.

I verified this with the data given but can not find the exact explanation why does it give one extra point after one relevant one?