Shooting method for ODE interpolation error

2 views (last 30 days)
Meva
Meva on 29 Oct 2014
Commented: Meva on 8 Nov 2014
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

Accepted Answer

Matt Tearle
Matt Tearle on 6 Nov 2014
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.
  1 Comment
Meva
Meva on 8 Nov 2014
Thank you so much Matt. I debugged and remove all unnecessary inputs, and that worked :)

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by