Inputs must be floats, namely single or double. (ODE45 function)

Hello,

I am getting this error when I use the ode45 solver:

" _Error using odearguments (line 113) Inputs must be floats, namely single or double.

Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in bvp4c (line 53) [phi, x] = ode45(@(phi,x) odefun(phi,x,y), phispan, xc)_"

My code:

y = - (5063413703368611*log((2*x)/(x + (x^2 - 711601842214333/9007199254740992)^(1/2))))/18014398509481984 - 
(5063413703368611*bessely(0, x))/18014398509481984  % Obtained by solving another equation with x as an independent variable

phispan = [beta:-h:phistar];  
[phi, x] = ode45(@(phi,x) odefun(phi,x,y), phispan, xc)
function dxdp = odefun(phi,x,y)
    dxdp = x*cos(phi)/(x*y-sin(phi))       
end

One thing to note. The code works fine when the y variable is not in the ODE. E.g:

phispan = [beta:-h:phistar];  
[phi, x] = ode45(@(phi,x) odefun(phi,x), phispan, xc)
function dxdp = odefun(phi,x)
    dxdp = x*cos(phi)/(x-sin(phi))       
end

Any help would be appreciated. Thank you!

Answer 1

[Walter Roberson]

You would not be able to sensibly define the y in that form without having defined

syms x

leading to a symbolic expression for y. So you are passing in a symbolic y to odefun, one involving the symbolic variable x.

And it appears that you are assuming that for any one call, the symbolic variable x in the symbolic y will be replaced with the numeric x that was passed in. It will not be -- they are different variables that live in different workspaces.

You should be defining

y = @(phi,x) - (5063413703368611*log((2*x)/(x + (x^2 - 711601842214333/9007199254740992)^(1/2))))/18014398509481984 - (5063413703368611*bessely(0, x))/18014398509481984  % Obtained by solving another equation with x as an independent variable

function dxdp = odefun(phi,x,y)
    dxdp = x*cos(phi)/(x*y(phi,x)-sin(phi));
end

Or alternately, assuming that the y has been calculated as a symbolic expression,

yf = matlabFunction(y);
[phi, x] = ode45(@(phi,x) odefun(phi,x,yf), phispan, xc)

with

function dxdp = odefun(phi,x,y)
    dxdp = x*cos(phi)/(x*y(x)-sin(phi));
end

Or at last resort, continue to call

[phi, x] = ode45(@(phi,x) odefun(phi,x,y), phispan, xc)

but with

function dxdp = odefun(phi,x,y)
    dxdp = x*cos(phi)/(x*double(subs(y))-sin(phi));
end

[Manojkumar Lokanathan]

Thank you very much for your help! I tried the 2nd and 3rd method, and both worked!