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

5 views (last 30 days)
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);
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)

Answers (1)

Walter Roberson
Walter Roberson on 27 Nov 2021
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.

Community Treasure Hunt

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

Start Hunting!