24小时热门版块排行榜

>论坛更新日志 (3520)
>虫友互识 (1027)
>硕博家园 (1021)
>导师招生 (181)
>休闲灌水 (127)
>考博 (119)
>博后之家 (104)
>海外博后 (100)
>文献求助 (96)
>学术会议 (73)
>论文投稿 (51)
>招聘信息布告栏 (47)
>催化 (35)
>考研 (13)
>基金申请 (12)
>电化学 (10)

返回列表

freebjx

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1284
散金: 50
帖子: 228
在线: 9.5小时
虫号: 612208
注册: 2008-09-25
性别: GG
专业: 物理学

[交流] 【求助】高斯积分语句问题已有2人参与

请问下面高斯积分语句有问题吗？
q(k)=intgauss(f,-inf,inf,8,0.1834346425,0.3626837834);

文献的结果应该是条曲线，我积出来的是一条直线

回复此楼

» 猜你喜欢

烦躁已经有5人回复

1楼 2010-04-23 11:29:13

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

xiegangmai

版主 (职业作家)

我没头衔

仿真EPI: 3
应助: 157 (高中生)
贵宾: 5.735
金币: 40609
散金: 2394
红花: 63
沙发: 4
帖子: 4046
在线: 1330.2小时
虫号: 653607
注册: 2008-11-13
性别: GG
专业: 机械动力学
管辖: 仿真模拟

freebjx(金币+1): 2010-04-23 15:15

你的f是什么啊？
用的哪个版本的MATLAB？
intgauss是你自己写的函数吗？我在2009a中没搜索到这个函数。

赞一下

回复此楼

明德厚学、求是创新

2楼2010-04-23 14:11:33

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

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

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

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

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖