필터 지우기
필터 지우기

Shooting method for ODE interpolation error

조회 수: 1 (최근 30일)
Meva
Meva 2014년 10월 29일
댓글: Meva 2014년 11월 8일
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

채택된 답변

Matt Tearle
Matt Tearle 2014년 11월 6일
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
Meva
Meva 2014년 11월 8일
Thank you so much Matt. I debugged and remove all unnecessary inputs, and that worked :)

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Matrix Indexing에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by