搜索内容 搜索用户

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

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类型，一会是数值，最好分开重命名呢 ...%本人计算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  E1 E3  P A1 A2 A3 e1 e3  p q1 q2 q3 q4 z K1 K2 K3 K4 K5 K;     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;    z=0:0.01:6.28;%相角大约为2*pi取两位小数     K1(z)=arrayfun(@(z)integral2(Mup,0.39,0.85,A3,A1,'Method','iterated','AbsTol',0,'RelTol',1e-6),z);     K2(z)=arrayfun(@(z)integral2(Mdown,3.9e-1,8.5e-1,A3,A1,'Method','iterated','AbsTol',0,'RelTol',1e-6),z);     K3(z)=arrayfun(@(z)integral2(Mup,3.9e-1,8.5e-1,A2,A3,'Method','iterated','AbsTol',0,'RelTol',1e-6),z);     K4(z)=arrayfun(@(z)integral2(Mdown,3.9e-1,8.5e-1,A2,A3,'Method','iterated','AbsTol',0,'RelTol',1e-6),z);     K(z)=(K1(z)./K2(z))-(K3(z)./K4(z));     plot(z,K(z));复制代码 改了下

15292166661 发表于 6 天前
 15292166661 发表于 2021-5-4 15:05 改了下 振幅表达式，也就是被积函数相当复杂，只能用数值积分了，新手一枚

15292166661 发表于 6 天前
 WarnerChang 发表于 2021-5-4 14:55 你这些变量名都有重复的！比如z一会是sym类型，一会是数值，最好分开重命名呢 ... 我不定义z为变量名，我后面怎么定义z的步长并画图呢
 您需要登录后才可以回帖 登录 | 注册 本版积分规则 回帖后跳转到最后一页