Intersection of a line with constant slope and a line with changing slopes

Hi, I have two different lines that represent different river profiles. The green line above sea-level (see jpg at the bottom of this post) represents the modern river channel whereas the yellow line represents what the channel may have looked like in the past. I'm trying to figure out where both lines cross. The code also has to be pretty general as I will be applying it to different river profiles, not just this specific one.

Here's specific informationa about the code:

idx = 2980; % calls upon river profile #2980 out of 10848
channel_height=0:20;
channel_length=channel_len(idx,:); % channel_len comes from dataset
plot(channel_length,channel_height), % plotting the upslope channel profile (green line)
hold on,
plot(-1000*shelf_len(idx,1:7),0:-20:-120), % plotting the shelf profile from dataset (red line, not important in this case), shelf_len comes from dataset
xlabel('distance (m)'),
ylabel('depth (m)')
yline(0,'--','sea-level')
channel_length_nanremoved=rmmissing(channel_length); % removing nans from channel data
positionofnan=find(isnan(channel_length)); % finding positions of nans in channel length
channel_height_lennanremoved=channel_height; % applying the length nan positions to channel height to remove their corresponding heights
channel_height_lennanremoved(positionofnan)=[];
shelf_len_nanremoved=-1000.*rmmissing(shelf_len(idx,1:7)); % removing nans from shelf to combine with channel

% Using a forloop to calculate slope between channel data points
n=2:length(channel_length_nanremoved);
for ii=2:length(channel_length_nanremoved)
    slope=(channel_height_lennanremoved(n)-channel_height_lennanremoved(n-1))./(channel_len_nanremoved(n)-channel_len_nanremoved(n-1));
end

% Generating the yellow line representing a previous position of the channel
m=mean(slope(end-9:end))
shelf_len_reversed=flip(shelf_len_nanremoved,2); % flipping shelf lengths as they start off more negative and get more positive
x=horzcat(shelf_len_reversed,channel_length_nanremoved); % combining shelf and channel lengths
x1=channel_length_nanremoved(end); % starting line from start of channel profile
position2=find(channel_length==channel_length_nanremoved(end));
y1=channel_height(position2);
y=(m*(x-x1))+y1; % determining y for given x
hplot=plot(x,y);

Answer 1

[dpb]

idx = 2980;                                 % calls upon river profile #2980 out of 10848
channel_height=0:20;
channel_length=channel_len(idx,:);          % channel_len comes from dataset
plot(channel_length,channel_height),        % plotting the upslope channel profile (green line)
hold on
plot(-1000*shelf_len(idx,1:7),0:-20:-120),  % plot shelf profile from dataset (red line, not important in this case), shelf_len comes from dataset
xlabel('distance (m)'),
ylabel('depth (m)')
yline(0,'--','sea-level')

ixData=~isnan(channel_length);              % logical vector nans in channel data
ch_length=channel_length(ixData);           % keep good data only -- length
ch_height=channel_height(ixData);           % and height
shelf_len_nanremoved=-1000.*rmmissing(shelf_len(idx,1:7)); % removing nans from shelf to combine with channel

% Calculate slope between channel data points (no loop needed)
slope=diff(ch_height)./diff(ch_length);

% Generating the yellow line representing a previous position of the channel
m=mean(slope(end-9:end));

OK, at the above is where I'd probably split ways (besides shortening variable names so much easier to scan the code, anyways)...

Where does the magic number "9" come from and why hardcoded into the code? Can't change anything easily that way.

Why use the average of the slopes instead of fitting a linear trend to the data? Should be nearly the same but a single outlier woudn't be as significant...

b=polyval(ch_length,ch_height,1);      % least square linear regression thru ch_length

Then to find the intersection use fsolve on the difference between the two lines.

ADDENDUM POST COMMENT:

OK, I'll concede using the mean of the means if what you really want...I'd suggest investigating the difference a little more though unless this is a respected technique in the application field.

% Generating the yellow line representing a previous position of the channel
nMeansUse=10;                               % at least make a variable can change
m=mean(slope(end-(nMeansUse_1):end));       % an average mean of means
x=[fliplr(ch_length),ch_length);            % combine shelf and channel lengths
x1=ch_length(end);                          % starting line from start of channel profile
ix2=find(channel_length==ch_length(end),1);
y1=channel_height(ix2);
y=(m*(x-x1))+y1;                            % determining y for given x
hplot=plot(x,y);

% find intersection of two...
fnUP=@(x) interp1(channel_length,channel_height,x);    % interpolated green line (upslope)
fnPR=@(x) polyval([m y1],x-x1);                        % evaluate yellow line (previous)
x0=fsolve(@(x) fnUP(x)-fnPR(x),x1);                    % try to find interesection of two

[Derrick Vaughn]

Hi dpb, thanks for your answer! I'm new to using MATLAB and I see now where some of my code could be condensed better. Thanks for those suggestions. For the number 9, I was wanting to draw the second, yellow river profile as representing what each channel would look like if I used the mean slope of the last 10 slopes of the current river profile. Sorry I didn't include those details in the question.

Is there any way I can use fsolve to look for the differences between the current river profile and what the river profile would look like using the mean of the last 10 slopes?

[Derrick Vaughn]

Hi dpb, everything works until I get to the last line of the code. I get this following error:

Error using griddedInterpolant

The coordinates of the input points must be finite values; Inf and NaN are not

permitted.

Error in interp1 (line 151)

F = griddedInterpolant(X,V,method);

Error in @(x)interp1(channel_length,channel_height,x)

Error in @(x)fnUP(x)-fnPR(x)

Error in fsolve (line 248)

fuser = feval(funfcn{3},x,varargin{:});

Caused by:

Failure in initial objective function evaluation. FSOLVE cannot continue.

[dpb]

Something's wrong in your input, then...the input for one or the other of the functions or X1 must have a NaN or Inf in it.

Not having specific data here, can't debug; all was air-code; you'll have to use your datasets to debug there...

ADDENDUM:

Reading the error text more thoroughly, the failure was actually in interp1 which is only used for fnUP the upper line. That tells you precisely that it appears that your x vector above includes either some missing values (NaN or Inf).

[Derrick Vaughn]

Hi dpb, first thanks again for all of the help. I really appreciate it. I fixed the NaN problem but the solution it gives me for x0 is for the last point in the river profile where they both start out (distance = 1.35 x 10^5). It doesn't give me the point where they cross (~1.0 x 10^5).

I also took out -x1 from the one of the lines in the code because I wasn't sure why it was there. It now reads: fnPR=@(x) polyval([m y],x);

This gave me an answer closer to the intersection, but not quite there. Could you elaborate a little why -x1 was included in that line of code?

[dpb]

y=(m*(x-x1))+y1;   

above; the origin was at x1,y1, not x===0. The anonymous function just incorporates that as well.

it's trying to find the intersection that would occur by extrapolating on out to the right by using the RH point as starting guess.

Need better initial guess within the region that brackets the intersection; evaluate the difference between the linear line and the observed points and find where the sign changes to bracket the point...and will tell you whether there is a solution or not if there isn't a sign change.

I did say it was air code didn't I??? :)

[Derrick Vaughn]

Hi dpb,

I've tried several things over the last few days but I think the closest I got to solving it was using the line of code below. However, it still hasn't gotten me the exact point where the lines cross and I'm not 100% convinced I went about it the right way based on your comment above. Can you provide a little more guidance, thanks!

x0=fsolve(@(x) fnUP(x)-fnPR(x),find(fnPR(y)<fnUP(y),1));

[dpb]

I don't have any data. Attach enough data/code to create the section trying to do the fit/intersection over, anyways

[Derrick Vaughn]

Hi dpb, I've updated the code (see below) to make the code length and variable names shorter and I've attached mat files (ch_len, ch_height, and sh_length) that contain data for the one river I am currently using in this example.

plot(ch_len,ch_height);
hold on
plot(shelf_len,0:-20:-120);
xlabel('distance (m)'),
ylabel('depth (m)')
yline(0,'--','sea-level');

slope=diff(ch_height)./diff(ch_len);
nMeansUse=10
m=mean(slope(end-(nMeansUse-1):end));      
x=[fliplr(sh_len),ch_len];            
x1=ch_length(end);                          
ix2=find(channel_len==ch_len(end),1);
y1=ch_height(ix2);
y=(m*(x-x1))+y1;                            
hplot=plot(x,y);
fnUP=@(x) interp1(ch_length,ch_height,x);
fnPR=@(x) polyval([m y],x-x1);
x0=fsolve(@(x) fnUP(x)-fnPR(x),find(fnPR(y)<fnUP(y),1))

[dpb]

Won't run standalone...

Unrecognized function or variable 'shelf_len'. 

> x1=ch_length(end);
The end operator must be used within an array index expression. 

Need something that will run after clear, I don't have anything that already exists in your environment.

[Derrick Vaughn]

Hi dpb, sorry about that. I thought I had fixed the entire code to reflect just the three mat files I attached above. Any part of the above code that says ch_length should just be ch_len and any part that says shelf_len should just be sh_len. Here's an updated code with those adjustments:

plot(ch_len,ch_height);
hold on
plot(sh_len,0:-20:-120);
xlabel('distance (m)'),
ylabel('depth (m)')
yline(0,'--','sea-level');
slope=diff(ch_height)./diff(ch_len);
nMeansUse=10
m=mean(slope(end-(nMeansUse-1):end));      
x=[fliplr(sh_len),ch_len];            
x1=ch_len(end);                          
ix2=find(ch_len==ch_len(end),1);
y1=ch_height(ix2);
y=(m*(x-x1))+y1;                            
hplot=plot(x,y);
fnUP=@(x) interp1(ch_len,ch_height,x);
fnPR=@(x) polyval([m y],x-x1);
x0=fsolve(@(x) fnUP(x)-fnPR(x),find(fnPR(y)<fnUP(y),1))

[dpb]

Oh. I was rushed for time so didn't catch just the variable name change was all...I'll try to look at over weekend; am preparing madly for local community college foundation annual meeting that we're having to do with virtual meeting this year instead so lots of new things to get done...

[dpb]

OK, looked at some this AM...two problems:

1. fnPR is missing the "1" on the "y1" intercept term in anonymous function definition; the bracketed term to polyval() are the linear fit coefficients slope and intercept(*); and
2. Starting point for fsolve was outside range of ch_len >= 0 so fnUP returns NaN when extrapolates; if use RH end, then the two points match and so that is a solution.

It seems that fsolve is much more sensitive to finding a solution of this system than I would have thought; probably the interpolation via interp1 is the problem in that it (fsolve) can't compute a smooth derivative and get past the inflection points...

>> [fnUP(ch_len).' fnPR(ch_len).' fnUP(ch_len).'-fnPR(ch_len).']
ans =
         0   -5.6492    5.6492
    1.0000   -3.1621    4.1621
    2.0000   -1.2879    3.2879
    3.0000   -0.0413    3.0413
    4.0000    1.4096    2.5904
    5.0000    2.8109    2.1891
    6.0000    4.5303    1.4697
    7.0000    5.3827    1.6173
    8.0000    6.2146    1.7854
    9.0000    6.4686    2.5314
   10.0000    8.6111    1.3889
   11.0000   10.3100    0.6900
   12.0000   12.0088   -0.0088
   13.0000   14.1163   -1.1163
   14.0000   15.0329   -1.0329
   15.0000   16.0546   -1.0546
   16.0000   16.8865   -0.8865
   17.0000   17.8936   -0.8936
   18.0000   18.4366   -0.4366
   19.0000   19.0000         0
>>

does show the two functionals (as corrected) do work as expected/desired; the third column is the difference which does show a zero crossing. I'd suggest in general going forward to do that calculation and find the location where the sign does change and use the value at that interval as the starting guess to fsolve. Or, of course, if there is no sign change, then the two don't intersect and you can flag that result as well.

With a better starting point, fsolve did find the interesection...

Above figure from

>> x0=fsolve(@(x) fnPR(x)-fnUP(x),1.E5)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>
x0 =
   9.6508e+04
>> hX=plot(x0,fnPR(x0),'Xr','markersize',10);

>> [fnPR(x0)-fnUP(x0)] % absolute error
ans =
    0.0019
>> 
>> ans/x0              % relative error
ans =
   1.9831e-08
>> 

(*) When you passed [m y] you built a polynomial for fnPR that was interpreted to be a 27th order polynomial of

fnPR(x)=m*(x-x1)^27 + y(1)*(x-x1)^26 + y(2)*(x-x1)^25 + ... + y(26)*(x-x1) + y(27)

Needless to say, the above function doesn't much represent the PR data... :)

MORAL: test the pieces of the puzzle independently; use the debugger if need to to poke around at what's going on. Too bad there isn't a way from polyval to be able to output the functional form of the polynomial, but it's just computational, not symbolic nor representational.

BTW, the dotted lines in the plot above are the result of evaluating the two functionals over ch_len; the first diagnostic I did to ensure they were indeed doing what thought/intended them to do. I had spotted the use of y for y1 initially; don't know if that had happened in earlier code you posted back and I missed it then or not...

[Derrick Vaughn]

Hi dpb, I got it to return the point today by finding the point when the sign of the difference changes. Thanks for all of your help!

[dpb]

Yes, that'll be the way to find a good starting point...I'm sorry but I don't have the time at the moment to try to find a more robust solution technique, but either end of the section with the sign change should be reliable I think if not quite as convenient as just throwing a dart somewhere in the general vicinity.