[未答复] 如何求解复杂嵌套函数定积分

[复制链接]
jookerer 发表于 2022-11-23 21:38:47
二重符号积分,被积函数(下划线部分)过于复杂而无法积出来,且对于积分变量tao11的积分上限仍为自身tao11,积分结果作为下一数值积分中的积分函数的一部分进一步积分。尝试使用matlabFunction转化为句柄函数但失败了。主要困扰是符号积分部分,想得到具体的积分数值结果,不知该如何做,求指点!!!:'(

syms z tao11 r1 s1 r s tao
%参数
c_l=268; c_s=230; rho_s=7180; rho_l=6980; epsilon=0.4;        
lamda_f=40; rho_f=1000; c_f=300;                             
c_p=0.5*(c_l+c_s);rho_p=0.5*(rho_s+rho_l);
c_pb=epsilon*c_p+(1-epsilon)*c_f; rho_pb=epsilon*rho_p+(1-epsilon)*rho_f;
g=9.8; hp=4.7*10^9; k=1; L=58.5*1000; q0=10; R=2;  rp=10*10^(-9); l=10;  Tin=500; Tout=498.27; Tm=505; Tini=20;
u=0.2; u_ref=0.5; omega=0.26*10^(-9); lamda_s=67; theta=0.2; tao1=0.14;
Pe_p = 16*R*l*u*u_ref*rho_f*c_f/lamda_f*(rp/(l*R))^2*(r/R);
  Pe_p = double(subs(Pe_p,r,1e-3));                          
  if Pe_p<0.67
     C0=3; C1=1.5;
  elseif (0.67<=Pe_p)&&(Pe_p<=250)
     C0=1.8; C1=0.18;
  else
     C0=3; C1=1/11;
  end
   lamda_b=lamda_f;
  lamda_e=lamda_b*(1+C0*epsilon*Pe_p^C1);
  alfa_e=lamda_e/(rho_pb*c_pb*u_ref*l);
  q = q0*L/lamda_e/Tin;

  c=c_l/c_s; beta=L/c_s/(Tm*(1-omega/rp)-((Tin-Tout)*s1+Tin)); bi=rp*hp/lamda_s; sigma=Tm*omega/rp/(Tm*(1-omega/rp)-((Tin-Tout)*s1+Tin));
  p1=0.5*(1+sigma/bi-beta/bi)-0.5*(beta-(c-1)*sigma)*sigma/bi^2+1/2/bi*(sigma-sigma/bi-1+bi)*(beta-(c-1)*sigma);
  q1=1/bi*(c-1)*sigma+2*(beta-(c-1)*sigma)*sigma/bi^2-(1+sigma/bi)-1/bi*(sigma-sigma/bi-1+bi)*(beta-(c-1)*sigma);
  c1=tao1-2.5*(beta-(c-1)*sigma)*sigma/bi^2+0.5*(1+sigma/bi-beta/bi)+1/2/bi*(sigma-sigma/bi-1+bi)*(beta-(c-1)*sigma);
  QL0=-Tm*(1-omega/rp)+(Tin-Tout)*s1+Tin+(-Tm*omega/rp-rp*hp/lamda_s*(Tm*(1-omega/rp)-(Tin-Tout)*s1-Tin))*((-q1+sqrt(q1^2-4*p1*(c1-lamda_s*L/(rho_s*c_s*rp^2*u_ref)*tao11)))/2/p1-1)+(-rp*hp/lamda_s*(Tm*(1-omega/rp)-(Tin-Tout)*s1-Tin)/2/(L/c_s-(c_l/c_s-1)*Tm*omega/rp)*(Tm*omega/rp+rp*hp/lamda_s*(Tm*(1-omega/rp)-(Tin-Tout)*s1-Tin))+rp*hp/lamda_s*(Tm*(1-omega/rp)-(Tin-Tout)*s1-Tin)*(1-rp*hp/lamda_s)-hp*Tm*omega/lamda_s)*((-q1+sqrt(q1^2-4*p1*(c1-lamda_s*L/(rho_s*c_s*rp^2*u_ref)*tao11)))/2/p1-1)^2;
  QL=-3*epsilon*hp/rp*QL0;
  F=l*QL/rho_pb/c_pb/u_ref/Tin+g*l*u/c_pb/Tin*sin(theta);    % F(s1,tao11)函数!!

%求根
       N = 5;                                    % 根的数量;
       incr = 3.14;                                  % 求根过程的增量;
    for i = 1 : N                          % 含有0阶贝塞尔函数的第2个及后面的零点; ksi>=1个根以上求解;
       mu(i) = fzero(@(z)z.*0.5*(besselj(-1,z)-besselj(1,z)),i*incr);
       ksi(i)= sqrt((mu(i)./R).^2+i^2*pi^2);
       f(i)=i*pi;      
      d(i) = integral(@(r1) 4./(R^2*J0(i).^2).*besselj(0,mu(i)/R.*r1),0,R).*2.*pi.*r1;
      Tb(i) = exp(-ksi(i).^2*k.*tao11)*d(i)*int(int(sin(f(i)*s1).*exp(ksi(i).^2*k*tao11).*F,s1,0,1),tao11,0,tao11);  
      T(i) = integral(matlabFunction(Tb(i)),0,2e-13)
      i=i+1;
   end

您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表