[已答复] 关于计算一个含参变量的二重积分的问题

[复制链接]
15292166661 发表于 7 天前
本帖最后由 15292166661 于 2021-5-3 17:11 编辑

%本人计算B介子CP破坏中,前后不对称的大小
% B+>pi+pi-pi+ 三体衰变计算代码
  %Wolfenstein parameters for CKM matrix elements from CKMfitter;
  lambda=0.2255;A=0.810;rhobar=0.145;etabar=0.343;
  rho=real((1-A^2.*lambda^4)^0.5.*(rhobar+1i.*etabar)/((1-lambda^2)^0.5.*(1-A^2.*lambda^4.*(rhobar+1i.*etabar))));
  eta=imag((1-A^2.*lambda^4)^0.5.*(rhobar+1i.*etabar)/((1-lambda^2)^0.5.*(1-A^2.*lambda^4.*(rhobar+1i.*etabar))));
   %CKM matrix elements from CKMfitter;
    Vud=1-0.5*lambda.^2;
    Vcd=-lambda+0.5*lambda^5*A.^2*(1-2*(rho+1i*eta));
    Vub=lambda^3.*A.*(rho-1i.*eta);
    Vts=-1*A.*lambda^2+A.*lambda^4.*(0.5-rho-1i.*eta);
    Vus=lambda;
    Vcb=A.*lambda^2-0.5.*lambda^8.*A.^3.*(rho.^2+eta.^2);
    Vcs=1-lambda.^2./2;
    Vtd=A.*lambda.^3.*(1-rho-1i.*eta);
    Vtb=1;
    %color number
    Nc=3;
      %Wilson coefficients, ref: 1303.的[12]
    c1=1.149;c2=-0.325;c3=0.021+0.0045*1i;c4=-0.045-0.0136*1i;c5=0.0134+0.0045*1i;
    c6=-0.056-0.0136*1i;
     %a_i's
    a1=c1+c2./Nc; a2=c2+c1./Nc;a3=c3+c4./Nc; a4=c4+c3./Nc;
    a5=c5+c6./Nc; a6=c6+c5./Nc;
    GF=1.16638e-5;
    mf0=0.5;mrho0=0.775;mB=5.270;mpi=0.139;mu=0.0023;md=0.0048;mb=0.0049;%Mev
    taorho=0.15;taof0=0.5;%rho和f0的衰变振幅,来源于张老师算稿
    F0=0.5;A0=0.28;F1=0.3;%A1形状因子对应的是B——》rho,F1形状因子对应的是B-》f0,来源于张老师的文章
    grhopipi=(48*pi*taorho/mrho0*(1-4*mpi^2/mrho0^2)^(3/2)).^0.5;%强相互作用耦合系数
    frho=0.2160;fpi=0.1304;%介子衰变常数
    gf0pipi=(48*pi*mf0*taof0/(1-4*mpi^2/mf0^2)^0.5).^0.5;

    syms S12 S13 ACP Mrho0 Mrho0bar Mf0 Mf0bar E1 E3  P x y z A1 A2 A3 e1 e3  p q1 q2 q3 q4 Mup Mdown;
    S12=0.1:0.1:25;deta=0.004*6.28:0.004*6.28:6.28;S13=0.1:0.1:25;


        E1=@(x)0.5.*sqrt(x);
    E3=@(x)(mB.^2-mpi.^2-x)./(2.*sqrt(x));
    P=@(x)(E1(x).^2-mpi.^2).*(E3(x).^2-mpi.^2);
    A1=@(x)2.*E1(x).*E3(x)+2*mpi^2+2.*sqrt(P(x));
    A2=@(x)2.*E1(x).*E3(x)+2*mpi^2-2.*sqrt(P(x));
    A3=@(x)(A1(x)+A2(x))./2;

      %Mf0=(2^0.5*GF*fpi*gf0pipi.*F0.*(mB.^2-mpi.^2)./(x-mf0^2+1i*mf0*taof0)).*(Vud.*conj(Vub)*a2-Vtd*conj(Vtb)*(a4+2*mpi^2*a6./((mu+md)*(mu+mb))));
     Mf0=@(x)(0.0002./(x-mf0^2+1i*mf0*taof0)).*( 0.3384 - 0.0548i);
   %Mrho0bar=(2^0.5*GF*grhopipi*mrho0./ (x-mrho0^2+1i*mrho0*taorho)).*(y-A3)*(Vub.*conj(Vud)*((-1/2^0.5)*a1*frho*F1+a2*fpi*A0)-Vtb*conj(Vtd)*(a4+2*(mpi^2)*a6*fpi*A0/((mu+md)*(mu+mb))));
   Mrho0bar=@(x,y)(6.22e-05./ (x-mrho0^2+1i*mrho0*taorho)).*(y-A3(x))*(   0.0104 + 0.0077i);
   %Mf0bar=(2^0.5*GF*fpi*gf0pipi.*F0.*(mB.^2-mpi.^2)./ (x-mrho0^2+1i*mrho0*taorho))*(Vub.*conj(Vud)*a2-Vtb*conj(Vtd)*(a4+2*(mpi^2)*a6/((mu+md)*(mu+mb))));
    Mf0bar=@(x)( 0.0002./ (x-mrho0^2+1i*mrho0*taorho))*( 0.2785 + 0.1997i);
    %Mrho=(2^0.5*GF*grhopipi*mrho0./(x-mf0^2+1i*mf0*taof0)).*(y-A3)*(Vud.*conj(Vub)*((-1/2^0.5)*a1*frho*F1+a2*fpi*A0)-Vtd*conj(Vtb)*(a4+2*mpi^2*a6*fpi*A0./((mu+md)*(mu+mb))))
    Mrho0=@(x,y)(6.22e-05./(x-mf0^2+1i*mf0*taof0)).*(y-A3(x))*( 0.0126 - 0.0022i);
    M=@(x,y)Mrho0(x,y).*exp(1*1i.*z)+Mf0(x);
    Mbar=@(x,y) Mrho0bar(x,y).*exp(1*1i.*z)+Mf0bar(x);
    Mup=@(x,y) abs(M(x,y)).^2-abs(Mbar(x,y)).^2;Mdown=@(x,y)abs(M(x,y)).^2+abs(Mbar(x,y)).^2;
    syms z;
   z=0:0.01:6.28;%相角大约为2*pi取两位小数
    K1(z)=arrayfun(@(z)integral2(Mup,0.39,0.85,A3,A1),z);
    K2(z)=arrayfun(@(z)integral2(Mdown,3.9e-1,8.5e-1,A3,A1),z);
    K3(z)=arrayfun(@(z)integral2(Mup,3.9e-1,8.5e-1,A2,A3),z);
    K4(z)=arrayfun(@(z)integral2(Mdown,3.9e-1,8.5e-1,A2,A3),z);
    K(z)=(K1(z)./K2(z))-(K3(z)./K4(z));

    plot(z,K(z));
不知为何,最后总是得出
错误使用 integral2Calc>integral2t/tensor (line 231)
输入函数必须返回 'double' 或 'single' 值。找到 'sym'。
二重积分
出错 integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

出错 integral2Calc (line 9)
    [q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

6 条回复


WarnerChang 发表于 6 天前
在命令窗口输入:
doc integral2查看帮助文档呢
integral2中函数定义时不需要函数名等定义为sym类型吧

15292166661 发表于 6 天前
WarnerChang 发表于 2021-5-4 14:34
在命令窗口输入:
doc integral2查看帮助文档呢
integral2中函数定义时不需要函数名等定义为sym类型吧 ...

您好,我把被积函数Mup,Mdown,移除了syms后面,但还是报错

WarnerChang 发表于 6 天前
15292166661 发表于 2021-5-4 14:48
您好,我把被积函数Mup,Mdown,移除了syms后面,但还是报错

你这些变量名都有重复的!比如z一会是sym类型,一会是数值,最好分开重命名呢

15292166661 发表于 6 天前
WarnerChang 发表于 2021-5-4 14:55
你这些变量名都有重复的!比如z一会是sym类型,一会是数值,最好分开重命名呢 ...
  1. %本人计算B介子CP破坏中,前后不对称的大小
  2. % B+>pi+pi-pi+ 三体衰变计算代码
  3.   %Wolfenstein parameters for CKM matrix elements from CKMfitter;
  4.   lambda=0.2255;A=0.810;rhobar=0.145;etabar=0.343;
  5.   rho=real((1-A^2.*lambda^4)^0.5.*(rhobar+1i.*etabar)/((1-lambda^2)^0.5.*(1-A^2.*lambda^4.*(rhobar+1i.*etabar))));
  6.   eta=imag((1-A^2.*lambda^4)^0.5.*(rhobar+1i.*etabar)/((1-lambda^2)^0.5.*(1-A^2.*lambda^4.*(rhobar+1i.*etabar))));
  7.    %CKM matrix elements from CKMfitter;
  8.     Vud=1-0.5*lambda.^2;
  9.     Vcd=-lambda+0.5*lambda^5*A.^2*(1-2*(rho+1i*eta));
  10.     Vub=lambda^3.*A.*(rho-1i.*eta);
  11.     Vts=-1*A.*lambda^2+A.*lambda^4.*(0.5-rho-1i.*eta);
  12.     Vus=lambda;
  13.     Vcb=A.*lambda^2-0.5.*lambda^8.*A.^3.*(rho.^2+eta.^2);
  14.     Vcs=1-lambda.^2./2;
  15.     Vtd=A.*lambda.^3.*(1-rho-1i.*eta);
  16.     Vtb=1;
  17.     %color number
  18.     Nc=3;
  19.       %Wilson coefficients, ref: 1303.的[12]
  20.     c1=1.149;c2=-0.325;c3=0.021+0.0045*1i;c4=-0.045-0.0136*1i;c5=0.0134+0.0045*1i;
  21.     c6=-0.056-0.0136*1i;
  22.      %a_i's
  23.     a1=c1+c2./Nc; a2=c2+c1./Nc;a3=c3+c4./Nc; a4=c4+c3./Nc;
  24.     a5=c5+c6./Nc; a6=c6+c5./Nc;
  25.     GF=1.16638e-5;
  26.     mf0=0.5;mrho0=0.775;mB=5.270;mpi=0.139;mu=0.0023;md=0.0048;mb=0.0049;%Mev
  27.     taorho=0.15;taof0=0.5;%rho和f0的衰变振幅,来源于张老师算稿
  28.     F0=0.5;A0=0.28;F1=0.3;%A1形状因子对应的是B——》rho,F1形状因子对应的是B-》f0,来源于张老师的文章
  29.     grhopipi=(48*pi*taorho/mrho0*(1-4*mpi^2/mrho0^2)^(3/2)).^0.5;%强相互作用耦合系数
  30.     frho=0.2160;fpi=0.1304;%介子衰变常数
  31.     gf0pipi=(48*pi*mf0*taof0/(1-4*mpi^2/mf0^2)^0.5).^0.5;

  32.     syms S12 S13  E1 E3  P A1 A2 A3 e1 e3  p q1 q2 q3 q4 z K1 K2 K3 K4 K5 K;
  33.     S12=0.1:0.1:25;deta=0.004*6.28:0.004*6.28:6.28;S13=0.1:0.1:25;


  34.         E1=@(x)0.5.*sqrt(x);
  35.     E3=@(x)(mB.^2-mpi.^2-x)./(2.*sqrt(x));
  36.     P=@(x)(E1(x).^2-mpi.^2).*(E3(x).^2-mpi.^2);
  37.     A1=@(x)2.*E1(x).*E3(x)+2*mpi^2+2.*sqrt(P(x));
  38.     A2=@(x)2.*E1(x).*E3(x)+2*mpi^2-2.*sqrt(P(x));
  39.     A3=@(x)(A1(x)+A2(x))./2;

  40.       %Mf0=(2^0.5*GF*fpi*gf0pipi.*F0.*(mB.^2-mpi.^2)./(x-mf0^2+1i*mf0*taof0)).*(Vud.*conj(Vub)*a2-Vtd*conj(Vtb)*(a4+2*mpi^2*a6./((mu+md)*(mu+mb))));
  41.      Mf0=@(x)(0.0002./(x-mf0^2+1i*mf0*taof0)).*( 0.3384 - 0.0548i);
  42.    %Mrho0bar=(2^0.5*GF*grhopipi*mrho0./ (x-mrho0^2+1i*mrho0*taorho)).*(y-A3)*(Vub.*conj(Vud)*((-1/2^0.5)*a1*frho*F1+a2*fpi*A0)-Vtb*conj(Vtd)*(a4+2*(mpi^2)*a6*fpi*A0/((mu+md)*(mu+mb))));
  43.    Mrho0bar=@(x,y)(6.22e-05./ (x-mrho0^2+1i*mrho0*taorho)).*(y-A3(x))*(   0.0104 + 0.0077i);
  44.    %Mf0bar=(2^0.5*GF*fpi*gf0pipi.*F0.*(mB.^2-mpi.^2)./ (x-mrho0^2+1i*mrho0*taorho))*(Vub.*conj(Vud)*a2-Vtb*conj(Vtd)*(a4+2*(mpi^2)*a6/((mu+md)*(mu+mb))));
  45.     Mf0bar=@(x)( 0.0002./ (x-mrho0^2+1i*mrho0*taorho))*( 0.2785 + 0.1997i);
  46.     %Mrho=(2^0.5*GF*grhopipi*mrho0./(x-mf0^2+1i*mf0*taof0)).*(y-A3)*(Vud.*conj(Vub)*((-1/2^0.5)*a1*frho*F1+a2*fpi*A0)-Vtd*conj(Vtb)*(a4+2*mpi^2*a6*fpi*A0./((mu+md)*(mu+mb))))
  47.     Mrho0=@(x,y)(6.22e-05./(x-mf0^2+1i*mf0*taof0)).*(y-A3(x))*( 0.0126 - 0.0022i);
  48.     M=@(x,y)Mrho0(x,y).*exp(1*1i.*z)+Mf0(x);
  49.     Mbar=@(x,y) Mrho0bar(x,y).*exp(1*1i.*z)+Mf0bar(x);
  50.     Mup=@(x,y) abs(M(x,y)).^2-abs(Mbar(x,y)).^2;Mdown=@(x,y)abs(M(x,y)).^2+abs(Mbar(x,y)).^2;

  51.    z=0:0.01:6.28;%相角大约为2*pi取两位小数
  52.     K1(z)=arrayfun(@(z)integral2(Mup,0.39,0.85,A3,A1,'Method','iterated','AbsTol',0,'RelTol',1e-6),z);
  53.     K2(z)=arrayfun(@(z)integral2(Mdown,3.9e-1,8.5e-1,A3,A1,'Method','iterated','AbsTol',0,'RelTol',1e-6),z);
  54.     K3(z)=arrayfun(@(z)integral2(Mup,3.9e-1,8.5e-1,A2,A3,'Method','iterated','AbsTol',0,'RelTol',1e-6),z);
  55.     K4(z)=arrayfun(@(z)integral2(Mdown,3.9e-1,8.5e-1,A2,A3,'Method','iterated','AbsTol',0,'RelTol',1e-6),z);
  56.     K(z)=(K1(z)./K2(z))-(K3(z)./K4(z));

  57.     plot(z,K(z));
复制代码

改了下

15292166661 发表于 6 天前

振幅表达式,也就是被积函数相当复杂,只能用数值积分了,新手一枚

15292166661 发表于 6 天前
WarnerChang 发表于 2021-5-4 14:55
你这些变量名都有重复的!比如z一会是sym类型,一会是数值,最好分开重命名呢 ...

我不定义z为变量名,我后面怎么定义z的步长并画图呢
您需要登录后才可以回帖 登录 | 注册

本版积分规则

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