Hello everyone. I am trying to calculate a double integral in matlab. However my program gave me "Conversion to double from function_handle is not possible" message. Please help me how to fix it. I am a beginner and I do not know how to do. Thank you.
clear all;
clc;
format long
syms x
syms y
syms ome1
syms ome2
syms ome12
syms z1
syms z2
syms z12
A = 1.85; %alpha
a = 2.656; % r_0
D = 0.269; %D_0
M_Fe = 56;
M_B = 27;
e = 2.71828;
h = 6.625*10^(-34);
h_bar = h/(2*pi);
k_B = 1.3824*10^(-23);
muy_Fe = M_Fe/(M_Fe+M_B);
muy_B = M_B/(M_Fe+M_B);
C_Fe = 0.8; %Ti le phan tram cua Fe
C_B = 0.2;
k_eff = (1+(5/3)*(muy_Fe^2+muy_B^2))*D*A^2;
k3 = -(1-muy_Fe^3+muy_B^3)*D*A^3;
k4 = (7/12+(203/324)*(muy_Fe^4+muy_B^4))*D*A^4;
M_FeB = (M_Fe+M_B)/2;
M = 2*C_B*M_FeB+(C_Fe-C_B)*M_Fe;
omega_D = 2*sqrt(k_eff/M);
theta_D = h_bar*omega_D/k_B;
sigma_01 = (3*h_bar*a/(2*pi))*(2*C_B*k3+(C_Fe-C_B)*k3)/(2*C_B*k_eff+(C_Fe-C_B)*k_eff);
sigma_02 = (h_bar*a/(2*pi))*(1/(2*C_B*k_eff+(C_Fe-C_B)*k_eff));
sigma_03 = (3*h_bar^2*a^2/(4*pi))*(2*C_B*k3+(C_Fe-C_B)*k3)/(2*C_Fe*k_eff+(C_Fe-C_B)*k_eff);
k0_eff = 0.25*M*omega_D^2;
p = 10;
Tmax = 1000;
dT = Tmax/p;
T = 0.000001;
C3 = zeros(1,p);
Tem = (0:dT:Tmax-dT)
for n = 1:p
beta = 1/(k_B*T);
ome1 = @(x) 2.*sqrt(k_eff/M).*abs(sin(a.*x./2));
ome2 = @(y) 2.*sqrt(k_eff/M).*abs(sin(a.*y./2));
ome12 = @(x,y) 2.*sqrt(k_eff/M).*abs(sin(a.*(x+y)./2));
z1 = exp(beta.*h_bar.*ome1(x));
z2 = exp(beta.*h_bar.*ome2(y));
z12 = exp(beta.*h_bar.*(ome1(x)+ome2(y)));
fun = @(x,y) (ome1(x).*ome2(y).*ome12(x,y)/(ome1(x)+ome2(y)+ome12(x,y)))*(1+(6.*(ome1(x)+ome2(y))./(ome1(x)+ome2(y)-ome12(x,y)))*((z1(ome1).*z2(ome2)-2.*z12(ome1,ome2))./((z1(ome1)-1).*(z2(ome2)-1).*(z12(ome1,ome2)-1))));
ymax = @(x) pi/a - x;
C3 = integral2(fun,0,1,0,ymax)
T = T+dT;
waitbar(n/p)
end
