magrog 发表于 2008-9-16 10:52

二重积分的问题

下面是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 = ;
    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 发表于 2008-9-16 11:34

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

magrog 发表于 2008-9-16 14:18

是将将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 发表于 2008-9-16 14:22

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

magrog 发表于 2008-9-16 14:37

回复 地板 sigma665 的帖子

=(((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 发表于 2008-9-17 10:50

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 发表于 2008-9-17 11:14

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

magrog 发表于 2008-9-17 13:57

当然是要把r=0.04去掉。

sigma665 发表于 2008-9-18 09:21

function result = myfun(x,r,h,A,d,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


clear all
clc
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,[],[],h,A1,d,r,fai); % 0~pi积分区间
temp= dblquad(@myfun,0,pi,0.04,0.05,[],[],h,A1,d,fai);
result = ;
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');

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

magrog 发表于 2008-9-18 10:26

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

sigma665 发表于 2008-9-18 10:34

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

个人习惯

magrog 发表于 2008-9-18 10:48

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

sigma665 发表于 2008-9-18 11:05

回复 12楼 magrog 的帖子

9楼可以运行

jojocleo 发表于 2008-9-18 12:19

双重积分慢有办法吗？

magrog 发表于 2008-9-18 18:11

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

[ 本帖最后由 magrog 于 2008-9-18 18:13 编辑 ]
页: [1] 2
查看完整版本: 二重积分的问题