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

clc; clear all; close all; 

h_bar = 1.054560652E-34/1.6E-19;              
k_b = 1.3806488E-23/1.6E-19;              


syms omega
d = 10E-9;   
a = 1;
tspan = [0.001,a];
f = linspace(0.01,0.1,5);
parfor y = 1:length(f)
    
    T_H = 450;            %T1 
    T_L = 300;            %T2
    
    Theta_omega_T_L = @(omega,kp) (h_bar*omega/2).*coth((h_bar*omega)/(2*k_b*T_L));
    Theta_omega_T_H = @(omega,kp) (h_bar*omega/2).*coth((h_bar*omega)/(2*k_b*T_H));
    THETA = @(omega,kp) Theta_omega_T_H(omega,kp) - Theta_omega_T_L(omega,kp);
    
    %..................Effective..............................%
    omegap = 2.7E14;
    eps2 = 1.77;
    eps_inf = 1;       
    gamma = 1E12; 
    k_om = 0.1;
    omegap_r = @(r)omegap.*(1-(k_om.*r./a));  
    eps = @(r) eps_inf - (omegap_r(r).^2./(omega.^2+(1i*omega.*gamma)));
    opts = odeset('RelTol',1.e-5,'AbsTol',1.e-8);
    [r,eps_bar] = ode45(@(r,eps_bar) (eps(r)-eps_bar).*(eps_bar+2.*eps(r))./(r.*eps(r)),tspan,1,opts);
    e_a = eps_bar(end);
    epsilon_e = @(omega) eps2+((3*eps2*f(y)).*((e_a(1)-eps2)./(e_a(1)+2*eps2)));
    
    kz = @(omega,kp) sqrt((omega/c).^2-kp.^2);
    k_1 = @(omega,kp) sqrt(epsilon_e(omega).*(omega/c).^2-kp.^2);                        
    
    
    r_1pp = @(omega,kp)(epsilon_e(omega).*kz(omega,kp)-k_1(omega,kp))./(epsilon_e(omega).*kz(omega,kp)+k_1(omega,kp));
    r_1ss = @(omega,kp)(kz(omega,kp)-k_1(omega,kp))./(kz(omega,kp)+k_1(omega,kp));
    
    
    xi_1 = @(omega,kp)((1 - abs(r_1ss(omega,kp)))^2);
    Xi_1 = @(omega,kp) xi_1(omega,kp).*kp.*THETA(omega,kp);
    ymin = @(omega) omega/3e8;
    XI_1 = integral2(Xi_1,1E13,1E15,1E-5,ymin,'RelTol',1E-3)
    
end

[Walter Roberson]

What is the difference between this and https://www.mathworks.com/matlabcentral/answers/1588714-i-am-having-error-using-ode45-with-inputs-must-be-floats-namely-single-or-double#comment_1856089 ?

Answer 1

[Walter Roberson]

syms omega

omega is symbolic.

    eps = @(r) eps_inf - (omegap_r(r).^2./(omega.^2+(1i*omega.*gamma)));

eps computes using omega, and omega is symbolic, so calling eps is going to return a symbolic value.

    [r,eps_bar] = ode45(@(r,eps_bar) (eps(r)-eps_bar).*(eps_bar+2.*eps(r))./(r.*eps(r)),tspan,1,opts);

the function there calls eps, but eps returns a symbolic value, so the function returns a symbolic value.

When you use ode45(), the function you give must not return a symbolic value.

There is no way to get around this problem when you use a numeric solver: you would have to use a symbolic solver and hope that the symbolic solver was able to produce a formula for you.

    e_a = eps_bar(end);
    epsilon_e = @(omega) eps2+((3*eps2*f(y)).*((e_a(1)-eps2)./(e_a(1)+2*eps2)));

Imagine for a moment that you were able to use dsolve() and so you got a symbolic solution for eps_bar, a solution in terms of omega. In such a case, e_a would be in terms of the symbolic variable omega. Inside epsilon_e, you use e_a so that would bring the symbolic omega into the expression. However, inside anonymous function bodies, named parameters, such as @(omega) are only substituted for direct references to the parameter inside the function. The symbolic omega that might hypothetically be output into e_a would not be the same as the name omega passed into the epsilon_e function.

Your epsilon_e does not refer to the passed-in omega, so if it were to work at all, it would return a constant value not affected by whatever was passed in.

If, hypothetically, you were able to get eps_bar out as an expression in terms of symbolic omega, then you would need to rewrite most of the rest of your code to deal symbolic expressions.

When I scan the rest of your code, it looks to me as if you are indeed hoping that ode45() can return an expression for eps_bar that is in terms of omega, and that you want to be able to do a 2D integral that includes numeric values for that omega.

If I am correct, then you are going to need to rewrite the ode into symbolic form and dsolve() it.

My tests suggest that doing that might indeed be possible.