Writing a linear interpolation function in MATLAB

Question

I'm trying to write a linear interpolation script in matlab. How do I fix that the vectors are of different length?

I've tried to wrap my head around this issue, but can't seem to get it. Should an if statement be included anywhere in the for loop?

z = linspace(0,1000,21)
vel = 1500*z^0.1;

% I want to interpolate vel between the 201 elements of depths. 
depths = [0:5:1000];
numel_depths = numel(depths);

for ii = 1:numel_depths   %Going through all indices of depths

    lower = find(z > depths(ii),1); % Finding x0 and x1.
    higher = find(z < depths(ii),1,'last'); 

    V2(ii) = vel(lower) + ((vel(higher) - vel(lower))/(z(higher)-  
        z(lower)))*(depths(ii)-z(lower)); % linear interpolation step
end

Now it returns an error saying that the different sides had different number of elements. Is there a way to fix this so that it works as the interp1 function already installed in MATLAB?

Are you using Octave? You cannot continue to the new line like that in MATLAB. Please provide the *exact wording* of the error message; not your interpretation of it. — Sardar Usama, Feb 01 '19 at 10:19
Are you trying to write this function as a learning exercise or because you want linear interpolation? If it's the latter, use `interp1`. — Wolfie, Feb 01 '19 at 10:32

Wolfie · Accepted Answer · 2019-02-01T16:04:40.880

There are several issues with your code:

You need the element-wise power operator to define vel, as instructed by the error message you must have seen:
```
vel = 1500*z.^0.1; % Use .^ rather than ^
```
There is a chance that lower or upper are empty, if no points are in the range you're testing. You need to skip over these cases using a test like:
```
~isempty( lower ) && ~isempty( higher )
```
You haven't used valid line-continuation. In MATLAB you cannot just break a line, you have to add ellipsis (...) to the end of the broken line.
You are using strict inequalities > and <, meaning you miss boundary points which you should be including in the interpolation (demo below).

You should pre-allocate arrays before for loops, so

v2 = NaN(size(depths)); % Define output V2 to be same size as input depths

lower is a built-in MATLAB function for making a string lower-case, avoid using it as a variable name.

Fixes to the above all come together like so:

z = linspace(0,1000,21);
vel = 1500*z.^0.1;      % Element-wise power operator
depths = 0:5:1000;
V2 = NaN(size(depths)); % Pre-allocate output array

for ii = 1:numel(depths);   
    % Finding x0 and x1, inclusive of boundaries (>= or <=).
    % Using 'x0' and 'x1' as var names to avoid shadowing the 'lower' function
    x0 = find(z >= depths(ii),1); 
    x1 = find(z <= depths(ii),1,'last'); 

    % Check if any points fell in this region
    if ~isempty( x0 ) && ~isempty( x1 )
        % Interpolation step, note the line continuation "..."
        V2(ii) = vel(x0) + ((vel(x1) - vel(x0))/(z(x1) - ...  
            z(x0)))*(depths(ii)-z(x0)); 
    end
end

We can validate this against the built-in interpolation function interp1

With your original strict inequalities < and >:

With non-strict inequalities as shown above <= and >=:

Writing a linear interpolation function in MATLAB

1 Answers1