Shooting method for ODE interpolation error

Hi,

I am using shooting method to solve second order ODE. The related part of main code is attached:

 function [f] = constriction_function(ii,dx,dt,c1,c2,H1D,AD,H1,A,fik,t,ub1,ub2,p2);

% d2T
% ---- = p, f(x=0)=f0, f(x=L)=f1
% dx2
% For the solution of the initial value problem
% the routine bvps is applied.
% z0: the initial steepness
f0 = 0; % boundary value on the left side
f1 = 0; % boundary value on the right side
% p2(nnn) % constant
L = 1; % length of the body
xs = 0; %coordinate of the initial point
f = f0; % initial condition
zstart=0.1;
deltaz=1;
Nz=1000;
zv(Nz+1)=zstart; 
z=zv(Nz+1);
[f,fvect,xvect]=co_constriction_function(ii,dx,dt,c1,c2,H1D,AD,H1,A,fik,t,ub1,ub2,p2);
fvegpont(Nz+1)=f;
for i=1:Nz
    zv(i)=zstart-(Nz+1-i)*deltaz; z=zv(i);
    [f,fvect,xvect]=co_constriction_function(ii,dx,dt,c1,c2,H1D,AD,H1,A,fik,t,ub1,ub2,p2);
    fvegpont(i)=f;
    zv(Nz+1+i)=zstart+i*deltaz; z=zv(Nz+1+i);
    [f,fvect,xvect]=co_constriction_function(ii,dx,dt,c1,c2,H1D,AD,H1,A,fik,t,ub1,ub2,p2);
    fvegpont(Nz+1+i)=f;
end
for i=1:2*Nz+1
    fvegpont(i); 
    zv(i);
end
% For finding the root we apply the inverse interpolation.
fint=fvegpont-f1; z=interp1(fint,zv,0);
fprintf('Steepness:  %fn',z)
[ffinal,fvect,xvect]=co_constriction_function(ii,dx,dt,c1,c2,H1D,AD,H1,A,fik,t,ub1,ub2,p2);

   end

There is a function the program uses but if I attached it, the question would be so long. I got the following error message:

Error using griddedInterpolant
Interpolation requires at least two sample points in each dimension.

Error in interp1 (line 186)
        F = griddedInterpolant(X,V,method);

Error in constriction_function (line 36)
fint=fvegpont-f1; z=interp1(fint,zv,0);

Actually my main program uses this function and another one. I searhed the error message but no luck. Any help please?

Best, Meva

Answer 1

[Matt Tearle]

Set a breakpoint on line 36 and run the code. When it stops at that line, see what you get from

unique(fvegpont)

I'm guessing you get a single value. So if you disp(fvegpont) or open up fvegpont in the Variable Editor, you'll see a vector of the same value ( f ).

As far as I can see, you make the exact same function call three times: [f,fvect,xvect]=co_constriction_function(ii,dx,dt,c1,c2,H1D,AD,H1,A,fik,t,ub1,ub2,p2); None of the inputs seems to change, so (a) this is a waste of computations, and (b) you have the same value for f (and fvect and xvect) each time. So when you do fvegpont(i)=f; and fvegpont(Nz+1+i)=f; in the loop, you're assigning the same value to every element.

You should also pay attention to the Code Analyzer warnings the Editor is giving you. There are a lot of variables that don't seem to be used at all. That's not a good sign.

[Meva]

Thank you so much Matt. I debugged and remove all unnecessary inputs, and that worked :)