Unknown variable as an index of a parameter inside of ode45

Hello!

I have the code shown below, in which an unknown variable Y(2) is used to compute an index in an array, passed as a parameter to ode45.

It looks like this is illegal, since I get an error:

Attempted to access cc(0); index must be a positive integer or logical.

Please provide some workaround. Thank you.

%%%%%%%%%%%%%%%%%%%%
function [s,y] = f(xo,zo, theta0,ss)
maxDepth = 5000;
zz = 0:1:maxDepth;
cc = 1520 + [zz*-.05];
cc(751:end) = cc(750) + (zz(751:end)-zz(750))*.014;
dcdz = diff(cc); % derivative of the sound speed w.r.t. z (depth)

%time interval and initial conditions
s_interval = [0,ss];
xio = cos(theta0)/cc(floor(zo)+1);
zetao = sin(theta0)/cc(floor(zo)+1);
init_cond = [xo,zo,xio,zetao]';
%solution
[s,y] = ode45(@(s,Y) rhs(s,Y,cc) , s_interval , init_cond);
%plot
plot(y(:,1), y(:,2), 'r');

function dYdt = rhs(s,Y,cc)
    Z = floor(Y(2))+1; Xi = Y(3); Zeta = Y(4);
dYdt = [ cc(Z)*Xi;
         cc(Z)*Zeta;
         0;
         -dcdz(Z)/cc(Z)^2];
end
end
%%%%%%%%%%%%%%%%%%%%

Answer 1

[Walter Roberson]

Your Y is your boundary conditions, and your boundary conditions are seldom going to be positive integers.

When I look at your code I suspect that what you are attempting to do is to interpolate into an array based upon the value of Y(2). If you were going to do that, then you should have those arrays and the marginal values and use interp1() .

However, unless you are pretty careful, using interpolation to calculate a coefficient for ode functions introduces discontinuities in one of the first two derivatives of your equations, which is a problem because the mathematics of RK45 requires continuity of the first two derivatives (it is sort of doing jacobian and hessian interpolation... sort of.)

It is possible to choose an interpolation function for interp1 that has the required continuity; in particular, cubic spline interpolation satisfies the condition mathematically.