| 查看: 805 | 回复: 5 | |||
freebjx金虫 (小有名气)
|
[交流]
【求助】高斯积分语句问题 已有2人参与
|
|
请问下面高斯积分语句有问题吗？ q(k)=intgauss(f,-inf,inf,8,0.1834346425,0.3626837834); 文献的结果应该是条曲线，我积出来的是一条直线 |
回复此楼» 猜你喜欢
小木虫没落了，除了祈祷帖子，几乎看不到有价值的帖子
已经有11人回复
-
会评
已经有3人回复
-
E口青C会评结束了吗？
已经有7人回复
-
生命口C13面上会评时间
已经有9人回复
-
内心匮乏
已经有13人回复
-
你可曾有过这种感受？
已经有6人回复
-
农业机械学报投稿周期
已经有7人回复
-
国自然面上D口祈祷
已经有4人回复
-
跨出版社商投稿
已经有4人回复
-
没收到消息，再次陪跑
已经有3人回复
2楼2010-04-23 14:11:33
freebjx
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1284
- 散金: 50
- 帖子: 228
- 在线: 9.5小时
- 虫号: 612208
- 注册: 2008-09-25
- 性别: GG
- 专业: 物理学
|
function I = IntGauss(f,a,b,n,AK,XK) if(n<5 && nargin == 4) AK = 0; XK = 0; else XK1=((b-a)/2)*XK+((a+b)/2); I=((b-a)/2)*sum(AK.*subs(sym(f),findsym(f),XK1)); end ta = (b-a)/2; tb = (a+b)/2; switch n case 0, I=2*ta*subs(sym(f),findsym(sym(f)),tb); case 1, I=ta*(subs(sym(f),findsym(sym(f)),ta*0.5773503+tb)+... subs(sym(f),findsym(sym(f)),-ta*0.5773503+tb)); case 2, I=ta*(0.55555556*subs(sym(f),findsym(sym(f)),ta*0.7745967+tb)+... 0.55555556*subs(sym(f),findsym(sym(f)),-ta*0.7745967+tb)+... 0.88888889*subs(sym(f),findsym(sym(f)),tb)); case 3, I=ta*(0.3478548*subs(sym(f),findsym(sym(f)),ta*0.8611363+tb)+... 0.3478548*subs(sym(f),findsym(sym(f)),-ta*0.8611363+tb)+... 0.6521452*subs(sym(f),findsym(sym(f)),ta*0.3398810+tb)... +0.6521452*subs(sym(f),findsym(sym(f)),-ta*0.3398810+tb)); case 4, I=ta*(0.2369269*subs(sym(f),findsym(sym(f)),ta*0.9061793+tb)+... 0.2369269*subs(sym(f),findsym(sym(f)),-ta*0.9061793+tb)+... 0.4786287*subs(sym(f),findsym(sym(f)),ta*0.5384693+tb)... +0.4786287*subs(sym(f),findsym(sym(f)),-ta*0.5384693+tb)+... 0.5688889*subs(sym(f),findsym(sym(f)),tb)); end |
3楼2010-04-23 15:17:06
freebjx
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1284
- 散金: 50
- 帖子: 228
- 在线: 9.5小时
- 虫号: 612208
- 注册: 2008-09-25
- 性别: GG
- 专业: 物理学
|
function [ output_args ] = Untitled10( input_args ) %UNTITLED10 Summary of this function goes here % Detailed explanation goes here clear clc r=1.0; c=3.0*10^8; h_b=(6.626196*10^-34)/(2*pi); k_b=1.3806505*10^-23; T=0.001; omegl=1.77*10^15; L=0.025; m=1.45*10^-10; kap=2.0*pi*215*10^3; omegm=2.0*pi*947*10^3; gamm=2.0*pi*140*1000; P=6.9*10^-3; N=sinh(r)*sinh(r); M=sinh(r)*cosh(r); omegc=2.0*10^15; gg=2.0*pi*2.7; omegm=2.0*pi*947000; E =5.5936e+26; nloop =101; deta0 = linspace(4.5*10^6,7.5*10^6,nloop); a= zeros(nloop,1); b= zeros(nloop,1); cs= zeros(nloop,1); q=zeros(nloop,1); deta=zeros(nloop,1); syms omeg for k=1:nloop a(k) =(1/6/E/gg^2*((omegm^2)*(-9*deta0(k)*omegm*kap^2-deta0(k)^3*omegm+27*E^2*gg^2+3*3^(1/2)*(omegm^2*kap^6+2*omegm^2*kap^4*deta0(k)^2+omegm^2*kap^2*deta0(k)^4-18*deta0(k)*omegm*kap^2*E^2*gg^2-2*deta0(k)^3*omegm*E^2*gg^2+27*E^4*gg^4)^(1/2)))^(1/3)-1/6*omegm^2*(3*kap^2-deta0(k)^2)./E/gg^2/(omegm^2*(-9*deta0(k)*omegm*kap^2-deta0(k)^3*omegm+27*E^2*gg^2+3*3^(1/2)*(omegm^2*kap^6+2*omegm^2*kap^4*deta0(k)^2+omegm^2*kap^2*deta0(k)^4-18*deta0(k)*omegm*kap^2*E^2*gg^2-2*deta0(k)^3*omegm*E^2*gg^2+27*E^4*gg^4)^(1/2)))^(1/3)+1/3/E/gg^2*deta0(k)*omegm)*kap; b(k) =(2*gg^2*E*(1/6/E/gg^2*(omegm^2*(-9*deta0(k)*omegm*kap^2-deta0(k)^3*omegm+27*E^2*gg^2+3*3^(1/2)*(omegm^2*kap^6+2*omegm^2*kap^4*deta0(k)^2+omegm^2*kap^2*deta0(k)^4-18*deta0(k)*omegm*kap^2*E^2*gg^2-2*deta0(k)^3*omegm*E^2*gg^2+27*E^4*gg^4)^(1/2)))^(1/3)-1/6*omegm^2*(3*kap^2-deta0(k)^2)/E/gg^2./(omegm^2*(-9*deta0(k)*omegm*kap^2-deta0(k)^3*omegm+27*E^2*gg^2+3*3^(1/2)*(omegm^2*kap^6+2*omegm^2*kap^4*deta0(k)^2+omegm^2*kap^2*deta0(k)^4-18*deta0(k)*omegm*kap^2*E^2*gg^2-2*deta0(k)^3*omegm*E^2*gg^2+27*E^4*gg^4)^(1/2)))^(1/3)+1/3/E/gg^2*deta0(k)*omegm)-deta0(k)*omegm)*(1/6/E/gg^2*(omegm^2*(-9*deta0(k)*omegm*kap^2-deta0(k)^3*omegm+27*E^2*gg^2+3*3^(1/2)*(omegm^2*kap^6+2*omegm^2*kap^4*deta0(k)^2+omegm^2*kap^2*deta0(k)^4-18*deta0(k)*omegm*kap^2*E^2*gg^2-2*deta0(k)^3*omegm*E^2*gg^2+27*E^4*gg^4)^(1/2)))^(1/3)-1/6*omegm^2*(3*kap^2-deta0(k)^2)/E/gg^2./(omegm^2*(-9*deta0(k)*omegm*kap^2-deta0(k)^3*omegm+27*E^2*gg^2+3*3^(1/2)*(omegm^2*kap^6+2*omegm^2*kap^4*deta0(k)^2+omegm^2*kap^2*deta0(k)^4-18*deta0(k)*omegm*kap^2*E^2*gg^2-2*deta0(k)^3*omegm*E^2*gg^2+27*E^4*gg^4)^(1/2)))^(1/3)+1/3/E/gg^2*deta0(k)*omegm)/omegm; cs(k)=a(k)+i*b(k); deta(k)=deta0(k)-2*gg^2*(cs(k)*cs(k)')/omegm; %-------------------------------------------- d1=-4.0*omegm*gg^2*deta(k)*(cs(k)*(cs(k)'))+(omegm^2-omeg^2-i*gamm*omeg)*((kap-i*omeg)^2+deta(k)^2); d2=-4.0*omegm*gg^2*deta(k)*(cs(k)*cs(k)')+(omegm^2-omeg^2+i*gamm*omeg)*((kap+i*omeg)^2+deta(k)^2); d3=-4.0*omegm*gg^2*deta(k)*(cs(k)*cs(k)')+(omegm^2- (2*omegm-omeg)^2 -i*gamm*(2*omegm-omeg))*((kap-i*(2*omegm-omeg))^2+deta(k)^2); d4=-4.0*omegm*gg^2*deta(k)*(cs(k)*cs(k)')+(omegm^2- (-2*omegm-omeg)^2 -i*gamm*(-2*omegm-omeg))*((kap-i*(-2*omegm-omeg))^2+deta(k)^2); A=(8*kap*gg^2*cs(k)*cs(k)'*((N+1)*( kap^2+(deta(k)+omeg)^2)+N* (kap^2+(deta(k)-omeg)^2))+(2*gamm*omeg/omegm)*((deta(k)^2+kap^2-omeg^2)^2+ 4*kap^2*omeg^2)*( 1+coth(h_b*omeg/(2*k_b*T))))/(d1*d2); B=8*kap*gg^2*(cs(k)')^2*M *(kap-i*(deta(k)+omeg))*(kap-i*(deta(k)+2*omegm-omeg))/(d1*d3); C=(8*kap*gg^2*(cs(k))^2*M')*(kap+i*(deta(k)-omeg))*( kap+i*(deta(k)+2*omegm+omeg))/(d1*d4); f=(omeg^2*A+omeg*(omeg-2*omegm)*B+omeg*(omeg+2*omegm)*C)/(2*pi); q(k)=intgauss(f,-inf,inf,8,0.1834346425,0.3626837834); if rem(k,10)==0, fprintf('%d ',k); end end plot(deta0/1e6,q); |
4楼2010-04-23 15:18:15
freebjx
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1284
- 散金: 50
- 帖子: 228
- 在线: 9.5小时
- 虫号: 612208
- 注册: 2008-09-25
- 性别: GG
- 专业: 物理学
5楼2010-04-23 16:32:02
freebjx
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1284
- 散金: 50
- 帖子: 228
- 在线: 9.5小时
- 虫号: 612208
- 注册: 2008-09-25
- 性别: GG
- 专业: 物理学
6楼2010-04-23 21:21:15
60