二重积分的问题

magrog

下面是myfun.m

function result = myfun(x,h,A,d,r,fai)
result =(((h+A*cos(x)*sin(fai)).^2).*(h^2+A^2+d^2+r^2+2*A*cos(x)*(h*sin(fai)-d*cos(fai))))./(((h^2+A^2+d^2+r^2+2*A*(h*sin(fai)-d*cos(fai))*cos(x))).^2-4*r^2*(d^2+A^2*(1-sin(fai)^2*cos(x).^2)-2*d*A*cos(fai)*cos(x))).^3/2;
end


下面是sjcll.m

i=1;
d=0.15;fai=pi/6;r=0.04;
for h=0.116:0.001:0.116
     result = [];
    for A1 = 0:0.0005:0.05
       temp= quad(@myfun,0,pi,[],0,h,A1,r,d,fai); % 0~pi积分区间
       result = [result,temp];
    end
    jie(i,:)=result;%把每次所得结果存成矩阵
     i=i+1;
end
A=0:0.0005:0.05;
for n=1:i-1
    plot(A,jie(n,:)/jie(n,1),'g.','Markersize',6)  
    hold on  
end
xlabel('A');ylabel('result');%标坐标
h1=0.05;
for j=1:i-1
    hh(j,1:length(A))=h1;               
    h1=h1+0.01;
end
figure
for m=1:i-1
    plot3(A,hh(m,:),jie(m,:),'r.','Markersize',6)
    hold on
end
xlabel('A');ylabel('h');zlabel('result');


上面的m文件对r的处理是赋值0.04，现在我想把原先的函数对x积分后再对r积分。积分范围从0.04到0.05.
试着将temp= quad(@myfun,0,pi,[],0,h,A1,r,d,fai); % 0~pi积分区间
改为： temp= dblquad(@myfun,0,pi,0.04，0.05，[],0,h,A1,d,fai); % 0~pi积分区间
可是运行一下，却发现行不通。
下面是报错的内容：
??? Error using ==> fcnchk
FUN must be a function, a valid string expression,
or an inline function object.
Error in ==> dblquad at 52
    quadf = fcnchk(quadf);
Error in ==> sjcll at 6
       temp= dblquad(@myfun,0,pi,0.04,0.05,[],0,h,A1,d,fai); % 0~pi积分区间

现在我是厚着脸皮来请教了。请多多包含。谢谢

sigma665

前一个能运行？
后一个不能运行？

magrog

是将将temp= quad(@myfun,0,pi,[],0,h,A1,r,d,fai); % 0~pi积分区间
改为： temp= dblquad(@myfun,0,pi,0.04，0.05，[],0,h,A1,d,fai); % 0~pi积分区间,
修改后的不能运行，不改的时候是正常的

sigma665

先说你方程中的未知数是什么
用word把方程打出来，然后再截图成jpg格式，发上来

magrog

=(((h+A*cos(x)*sin(fai)).^2).*(h^2+A^2+d^2+r^2+2*A*cos(x)*(h*sin(fai)-d*cos(fai))))./(((h^2+A^2+d^2+r^2+2*A*(h*sin(fai)-d*cos(fai))*cos(x))).^2-4*r^2*(d^2+A^2*(1-sin(fai)^2*cos(x).^2)-2*d*A*cos(fai)*cos(x))).^3/2;
上面的是函数。里面的h,A,d,fai,为参数，通过赋值进行积分。现在r,x为积分变量，先对x积分，然后对r积分。
没发现从哪里贴图。

[ 本帖最后由 magrog 于 2008-9-16 14:49 编辑 ]

magrog

temp= quad(@myfun,0,pi,[],0,h,A1,r,d,fai); % 0~pi积分区间
关键还是在这里，改为二次积分后temp= dblquad(@myfun,0,pi,[], 【】,h,A1,r,d,fai); % 0~pi积分区间
但是这样改，还是有问题。请高手们看看，能否解决啊？

sigma665

r=0.04前面已经定义了
直接改不行。。。

magrog

当然是要把r=0.04去掉。

sigma665

1. function result = myfun(x,r,h,A,d,fai)
   
2. result =(((h+A*cos(x)*sin(fai)).^2).*(h^2+A^2+d^2+r^2+2*A*cos(x)*(h*sin(fai)-d*cos(fai))))./(((h^2+A^2+d^2+r^2+2*A*(h*sin(fai)-d*cos(fai))*cos(x))).^2-4*r^2*(d^2+A^2*(1-sin(fai)^2*cos(x).^2)-2*d*A*cos(fai)*cos(x))).^3/2;
   
3. end
   

复制代码

1. clear all
   
2. clc
   
3. i=1;
   
4. d=0.15;fai=pi/6;%r=0.04;
   
5. for h=0.116:0.001:0.116
   
6. result = [];
   
7. for A1 = 0:0.0005:0.05
   
8. %temp= quad(@myfun,0,pi,[],[],h,A1,d,r,fai); % 0~pi积分区间
   
9. temp= dblquad(@myfun,0,pi,0.04,0.05,[],[],h,A1,d,fai);
   
10. result = [result,temp];
    
11. end
    
12. jie(i,:)=result;%把每次所得结果存成矩阵
    
13. i=i+1;
    
14. end
    
15. A=0:0.0005:0.05;
    
16. for n=1:i-1
    
17. plot(A,jie(n,:)/jie(n,1),'g.','Markersize',6)
    
18. hold on
    
19. end
    
20. xlabel('A');ylabel('result');%标坐标
    
21. h1=0.05;
    
22. for j=1:i-1
    
23. hh(j,1:length(A))=h1;
    
24. h1=h1+0.01;
    
25. end
    
26. figure
    
27. for m=1:i-1
    
28. plot3(A,hh(m,:),jie(m,:),'r.','Markersize',6)
    
29. hold on
    
30. end
    
31. xlabel('A');ylabel('h');zlabel('result');

复制代码

请注意参数的顺序
PS:双重积分运行时间较长

magrog

谢谢。昨天我把所有的A换成了a（偶尔看到大写字母表示矩阵，不知道是不是这样），在一定范围内竟然能够运算。因为运算时间很长，我就把A的范围缩小，竟然报错。真的是很奇怪。
谢谢楼上小西主任！

sigma665

昨天我把所有的A换成了a

个人习惯

magrog

问一下小西，上面的代码你运行过吗？我昨天运行了一下，好像是会报错的

sigma665

9楼可以运行

jojocleo

双重积分慢有办法吗？

magrog

d和h取相同数值时，运行出错，或在某些范围内也是出错。感觉很奇怪。是代码的问题还是函数在积分的时候就是有奇异点之类的东西呢？
Warning: Maximum function count exceeded; singularity likely.

[ 本帖最后由 magrog 于 2008-9-18 18:13 编辑 ]