-2

I have astraightline with the following data

xinter=[1.13 1.36 1.62 1.81 2.00 2.30 2.61 2.83 3.05 3.39]
yinter=[0.10 0.25 0.40 0.50 0.60 0.75 0.90 1.00 1.10 1.25]

and I want to find the intersection with a result of an interpolated data of a curve such as below

      a50=  [0.77 0.73 0.77 0.85 0.91 0.97 1.05 1.23 1.43 1.53 1.62 1.71 1.89 2.12 2.42];
a25=  [0.51 0.60 0.70 0.80 0.85 0.90 0.96 1.09 1.23 1.30 1.36 1.41 1.53 1.67];
vel25=[0.43 0.35 0.30 0.27 0.25 0.24 0.22 0.21 0.22 0.24 0.25 0.27 0.30 0.35];
vel50=[0.68 0.57 0.49 0.43 0.40 0.38 0.36 0.34 0.36 0.38 0.40 0.43 0.49 0.57 0.68 ];


% back up original data, just for final plot
bkp_a50 = a50 ; bkp_vel50 = vel50 ;

% make second x vector monotonic
istart = find( diff(a50)>0 , 1 , 'first') ;
a50(1:istart-1) = [] ;
vel50(1:istart-1) = [] ;

% prepare a 3rd dimension vector (from 25 to 50)
T = [repmat(25,size(a25)) ; repmat(40,size(a50)) ] ;
% merge all observations together
A = [  a25 ;   a50] ;
V = [vel25 ; vel50] ;

% find the minimum domain on which data can be interpolated
% (anything outside of that will return NaN)
Astart = max( [min(a25) min(a50)] ) ;
Astop  = min( [max(a25) max(a50)] ) ;

% use the function 'griddata'
[TI,AI] = meshgrid( 25:40 , linspace(Astart,Astop,10)  ) ; 
VI = griddata(T,A,V,TI,AI) ;

% plot all the intermediate curves
%plot(AI,VI)
hold on
% the original curves
%plot(a25,vel25,'--k','linewidth',2)
%plot(bkp_a50,bkp_vel50,'--k','linewidth',2)
% Highlight the curve at T = 30 ;
c30 = find( TI(1,:) == 40 ) ;
plot(AI(:,c30),VI(:,c30),'--r','linewidth',2)

xinter=[1.13 1.36 1.62 1.81 2.00 2.30 2.61 2.83 3.05 3.39]
yinter=[0.10 0.25 0.40 0.50 0.60 0.75 0.90 1.00 1.10 1.25]


x1inter=(AI(:,c30))';
y1inter=(VI(:,c30))';


yy2 = interp1(xinter, yinter, x1inter,'spline')


plot(xinter,yinter, '--k','linewidth',2)

 idx = find((y1inter - yy2) < eps, 1); %// Index of coordinate in array
 px = x1inter(idx)
py = y1inter(idx)
plot(px, py, 'ro', 'MarkerSize', 18)

But there is an error in the result when I modify x1inter

user183060
  • 13
  • 5

1 Answers1

-1

You can use piecewise polynomial curvefitting and the fzero function to find the intersection point:

pp1 = pchip(xinter,yinter);         % Curve 1
pp2 = pchip(AI(:,c30),VI(:,c30));   % Curve 2
fun = @(x) ppval(pp1,x) - ppval(pp2,x); % Curve to evaluate
xzero = fzero(fun,mean(xinter)) % intersection x value
yzero = ppval(pp1,xzero)
plot(xzero, yzero, 'bo', 'MarkerSize', 18)
Adriaan
  • 17,741
  • 7
  • 42
  • 75
M_Tornack
  • 114
  • 6
  • the result is not excatly producing the result but its very close it should be xzero =1.36, yzero=.25 – user183060 Dec 16 '16 at 12:01
  • I also zoomed in your plot and the marker is not really in the intersection point. Thanks – user183060 Dec 16 '16 at 12:02
  • The marker is not on the intersection point, because your plot uses linear interpolation. The calculated intersection point however is calculated with the Piecewise Cubic Hermite Interpolating Polynomial. Your intersection point (1.36 / 0.25) is given somewhere and you have to prove it? – M_Tornack Dec 16 '16 at 13:18
  • its working only the lower point of the interpolated data using c30 = find( TI(1,:) == 25) ; but if I use c30 = find( TI(1,:) == 40 ) ; it doesn't produce the exact intersection. It should be the xzero and yzero =(1.36, .25) for c30 = find( TI(1,:) == 25) and (1.62,0.40) as in the data – user183060 Dec 16 '16 at 16:55
  • Hello, I think it is a question what is for you the "exact intersection". Do you want to have the point on a discrete datapoint or a interpolation between 2 datapoints? As far as I can see is the point [1.36 | 0.25] a discrete point on the black curve. But this point is not the intersection point (but it is the closest). I chose interpolation to calculate the intersection point and get [1.38 | 0.26] for c30 = find( TI(1,:) == 25) and [1.64 | 0.41] for c30 = find( TI(1,:) == 40 ). – M_Tornack Dec 19 '16 at 08:06
  • Yes I want the exact interpolation for max value c30 = find( TI(1,:) == 40 ) which should be the px and py =(1.36, .25), and for the minimum value c30 = find( TI(1,:) == 25) is (1.62,0.40). And for c30 = find( TI(1,:) == 30) should be equal to the interpolated intersection of curve 1 and straight line – user183060 Dec 19 '16 at 09:12
  • please refer to the background about the in initial part of the code http://stackoverflow.com/questions/40989753/plot-interpolate-a-curve-between-2-different-types-of-curves-in-matlab – user183060 Dec 19 '16 at 09:13