Input function must return 'double' or 'single' values. Found 'sym'.

1 次查看（过去 30 天）

mohamed elkasemy 2021-10-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1563676-input-function-must-return-double-or-single-values-found-sym

clc;
clear all; 
close all;
C = 3*10^10;
KT = 0.086173324;  % MeV
M0 = 1.0078;       % u
M1 = 1.0078;       % u
M2 = 2.014102;     % u
M3 = 1.0078;       % u
M4 = 2.014102;     % u
M5 = 2.014102;     % u
M6 = 3.016030;     % u
M7 = 2.014102;     % u
M8 = 3.016030;     % u
M9 = 3.016030;     % u
M10 = 3.016030;    % u
M11 = 4.002603;    % u
Z0 = 1;  
Z1 = 1;  
Z2 = 1;  
Z3 = 1;
Z4 = 1;  
Z5 = 1;  
Z6 = 2;  
Z7 = 1;
Z8 = 2;  
Z9 = 2;
Z10 = 2;
Z11 = 2;
Eat = 0.98951013;
m01 = (M0*M1)/(M0+M1);
m23 = (M2*M3)/(M2+M3);
m45 = (M4*M5)/(M4+M5);
m67 = (M6*M7)/(M6+M7);
m89 = (M8*M9)/(M8+M9);
m1011 = (M10*M11)/(M10+M11);
cont = (931.494)/(3*10^10)^2;
cont1 = 10^-24;
syms y(t) 
sy0 =5; 
Q1 = @(t)(((0.6*10^-9)*((2*1.202)./(pi).^2)*(((10^3).*8065.73).^3))/(2*(1+exp((-1.293)./(t)))))/(6.023*10^23);
Q2 = @(t)(sqrt(8/(pi*m01*cont))).*((6.023*10^23)/(t)^(3/2)).*(integral (@(x) (exp(-Eat*Z0*Z1*sqrt((m01)./(x)))).*((3.94*10^-25) + (4.61*10^-24)*(x) + (2.96*10^-23)*((x).^2)).*cont1.*(exp(-(x)/(t))), 0,inf));
Q3 = @(t)(((0.6*10^-9).*((2*1.202)./(pi).^2)*(((10^3).*8065.73).^3))/(1*(1+exp((-1.293)./(t)))))/(6.023*10^23);
Q4 = @(t)(sqrt(8/(pi*m23*cont))).*((6.023*10^23)/(t)^(3/2)).*(integral (@(x) (exp(-Eat*Z2*Z3*sqrt((m23)./(x)))).*((0.20*10^-6) + (5.60*10^-6)*x + (3.10*10^-6)*(x.^2)).*cont1.*(exp(-(x)/(t))), 0,inf));
WV2   =@(t)(3*10^10)*sqrt(1-(1/((t/0.511)+1))^2);
WNep2 =@(t)((((3*1.202)./(2*pi.^2))*(((10^6)*sy0.*8065.73).^3)))*(exp((-0.511)./(t))); 
WNp2  =@(t)((0.6*10^-9)*((2*1.202)./(pi).^2)*(((10^6)*sy0.*8065.73).^3))/(1+exp((-1.293)./(t)));   
WNen2 =@(t)((((3*1.202)./(2*pi.^2))*(((10^6)*t.*8065.73).^3)))+ WNp2(t);
WNn2  =@(t) WNp2(t)*(exp((-1.293)./(t)));
Wtn2  = 1000;
Wconst2 = ((1.239*10^-10)^3); 
Wb2 =((7.56*10^-16)*((1.239*10^-6)^3))/((1.6*10^-19)*((8.617*10^-5)^4));
WF1 =@(t)(((1.166*10^-11)/(8065.73*10^6))^2)*((t)^2);
WF2 =@(t)(((1.166*10^-11)/(8065.73*10^6))^2)*((t)^2)*(exp(-(0.783)./t));    
WUU =@(t) (-(((Wconst2).*((WNep2(t).*WNn2(t).*WF1(t).*WV2(t))+(WNen2(t).*WNp2(t).*WF2(t).*WV2(t))+(WNep2(t).*WNen2(t).*WF1(t).*WV2(t))+(WNn2(t)./Wtn2)))/(4*Wb2.*((t)^2))));
ode(t) = diff(y(t),t) == ((Q1(t)*Q2(t))-(Q3(t)*Q4(t)*y(t)))/WUU(t);
cond = y(0) == 10^-5;
ySol(t) = dsolve(ode,cond);
t = linspace(0.01,10,100000);
loglog(t,ySol(t),'k-','linewidth',2)