请教多重积分问题 谢谢大家

weilinhy

我用matlab做上面的积分
1）符号积分
syms T H theta x y z;
int(int(int(exp(-theta*x-theta*abs(y-z))*(abs(x-y))^(2*H-2)*z^(2*H-2),x,0,T),y,0,T),z,0,T)
结果：Warning: Explicit integral could not be found.
2）数值积分
theta=1; T=10; H=0.78;
triplequad(inline('exp(-theta*x-theta*abs(y-z))*(abs(x-y))^(2*H-2)*z^(2*H-2)','x','y','z'),  0, T, 0, T, 0, T,1e-7,@quadl)
运行结果：
Error using ==> inlineeval
Error in inline expression ==> exp(-theta*x-theta*abs(y-z))*(abs(x-y))^(2*H-2)*z^(2*H-2)
??? Error using ==> eval
Undefined function or variable 'theta'.
Error in ==> inline.subsref at 25
    INLINE_OUT_ = inlineeval(INLINE_INPUTS_, INLINE_OBJ_.inputExpr, INLINE_OBJ_.expr);
我想请问
1）这个积分matlab可以符号积分表示出来么？
2）如果不行，数值积分求解应该怎么求呢！
谢谢大家

maigicku

关于出错信息的解释与改正请参考下帖：
常见的程序出错问题整理

1) 由上帖的解释可以看出，积分可能没有显式解，因而不适合用符号积分，改为数值积分;
2) 数值积分可以用匿名函数求解

weilinhy

匿名函数怎么弄呢？拜谢阿！刚接触matlab  很多不懂  谢谢了

maigicku

举个例子：
clear all
kk=linspace(0,5);
y=zeros(size(kk));
ff=@(k) ['sin(',num2str(k),'*x).*x.^2'];
f=@(k) quadl(ff(k),1,5);
for ii=1:length(kk)
y(ii)=f(kk(ii));
end
plot(kk,y)

上面例子中，@后面的(k)表示输入参数，ff就被表示成关于k的函数，同样，f也被表示为关于K的函数，这样就可以避免使用符号积分或内嵌函数。。如果不懂的话可以去搜一下，或是找本Matlab7的书。

weilinhy

谢谢maigicku
请问这个运行
theta=1; T=10; H=0.78;
ff=@(x,y,z) exp(-theta*x-theta*abs(y-z))*(abs(x-y))^(2*H-2)*z^(2*H-2);
f=@(k) quadl(ff(x,y,z),0,T,0,T,0,T);
结果为：Error: "x" is not the name of a function nor a variable,
but is used in an anonymous function either at the prompt or in the argument of EVAL.
帮忙看看 拜谢了

maigicku

quadl使用错误。。关于多重积分请参考下帖：
http://forum.vibunion.com/forum/vi ... hlight=%B6%FE%D6%D8
http://forum.vibunion.com/forum/viewthread.php?tid=83426&highlight=%B6%FE%D6%D8
f=quadl(@(z) quadl(@(y) quadl(@(x) exp(-theta*x-theta*abs(y-z))*(abs(x-y))^(2*H-2)*z^(2*H-2),0,T),0,T),0,T)

[ 本帖最后由 maigicku 于 2010-4-16 08:38 编辑 ]

weilinhy

input:
theta=1; T=10; H=0.78;
f=quadl(@(z) quadl(@(y) quadl(@(x) exp(-theta.*x-theta.*abs(y-z)).*(abs(x-y))^(2*H-2).*z^(2*H-2),0,T),0,T),0,T)

result:
??? Error using ==> mpower
Matrix must be square.

Error in ==> @(x)exp(-theta.*x-theta.*abs(y-z)).*(abs(x-y))^(2*H-2).*z^(2*H-2)

Error in ==> quadl at 70
y = feval(f,x,varargin{:}); y = y(:).';

Error in ==> @(y)quadl(@(x)exp(-theta.*x-theta.*abs(y-z)).*(abs(x-y))^(2*H-2).*z^(2*H-2),0,T)

Error in ==> quadl at 70
y = feval(f,x,varargin{:}); y = y(:).';

Error in ==> @(z)quadl(@(y)quadl(@(x)exp(-theta.*x-theta.*abs(y-z)).*(abs(x-y))^(2*H-2).*z^(2*H-2),0,T),0,T)

Error in ==> quadl at 70
y = feval(f,x,varargin{:}); y = y(:).';

[ 本帖最后由 ChaChing 于 2010-6-22 23:36 编辑 ]

rocwoods

楼主，原被积函数含有绝对值，需要按y和z，x和y的关系拆成四个积分的和，解决代码如下（需要R2009a以上的版本）：

1. 
   
2. clear all
   
3. theta=1; T=10; H=0.78;
   
4. ff1=@(x) @(y,z) exp(-theta*x-theta*(y-z)).*((x-y)).^(2*H-2).*z.^(2*H-2);
   
5. ff2=@(x) @(y,z) exp(-theta*x-theta*(y-z)).*((y-x)).^(2*H-2).*z.^(2*H-2);
   
6. ff3=@(x) @(y,z) exp(-theta*x-theta*(z-y)).*((x-y)).^(2*H-2).*z.^(2*H-2);
   
7. ff4=@(x) @(y,z) exp(-theta*x-theta*(z-y)).*((y-x)).^(2*H-2).*z.^(2*H-2);
   
8. Int1 = quadgk(@(x) arrayfun(@(x) quad2d(ff1(x),0,x,0,@(y)y),x),0,T)
   
9. Int2 = quadgk(@(x) arrayfun(@(x) quad2d(ff2(x),x,T,0,@(y)y),x),0,T)
   
10. Int3 = quadgk(@(x) arrayfun(@(x) quad2d(ff3(x),0,x,@(y)y,T),x),0,T)
    
11. Int4 = quadgk(@(x) arrayfun(@(x) quad2d(ff4(x),x,T,@(y)y,T),x),0,T)
    
12. INT = Int1+Int2+Int3+Int4
    
13. 
    
14. Int1 =
    
15. 
    
16.     1.2617
    
17. 
    
18. 
    
19. Int2 =
    
20. 
    
21.     4.1067
    
22. 
    
23. 
    
24. Int3 =
    
25. 
    
26.     1.4130
    
27. 
    
28. 
    
29. Int4 =
    
30. 
    
31.     3.2933
    
32. 
    
33. 
    
34. INT =
    
35. 
    
36.    10.0747
    

复制代码

关于上述代码的解释，我即将出版的书里有较多的篇幅进行介绍。这在现有的MATLAB参考书里很难见到。

[ 本帖最后由 rocwoods 于 2010-4-16 11:23 编辑 ]

weilinhy

多谢 多谢 rocwoods

maigicku

rocwoods 大神出现~~
LZ不好意思，我只不过给了个形式，并没有进行实际调试。。
还有rocwoods 大神：
ff1=@(x) @(y,z) exp(-theta*x-theta*(y-z)).*((x-y)).^(2*H-2).*z.^(2*H-2);
这种命名是不是只有R2009a才可以这样？
我还有个问题，能不能帮忙解决一下：
http://forum.vibunion.com/forum/thread-91762-1-2.html

rocwoods

呵呵，别叫大神，大家都是MATLAB学习者：）
那种用法不是R2009a特有，这是双重匿名函数的用法。7.0以上的版本都支持。主要是quad2d函数是2009a才有的，上面用到这个函数了。不用这个函数也可以，就是稍微麻烦些，要用两次arrayfun，书里面有介绍。
拟合的问题可以问下xiezhh，谢老师很在行的。呵呵

weilinhy

给您添麻烦了，谢谢您的回答

weilinhy

求解上面积分问题, 参考rocwoods给的方法
clear all; theta=1; T=10; H=0.78;
ff1=@(s2,u2)@( s1,u1)exp(-theta*(s2-u2)-theta*(s1-u1)).*(s2-s1).^(2*H-2).*(u2-u1).^(2*H-2);
ff2=@(s2,u2)@( s1,u1)exp(-theta*(s2-u2)-theta*(s1-u1)).*(s2-s1).^(2*H-2).*(u1-u2).^(2*H-2);

ff3=@(s2,u2)@( s1,u1)exp(-theta*(s2-u2)-theta*(s1-u1)).*(s1-s2).^(2*H-2).*(u2-u1).^(2*H-2);
ff4=@(s2,u2)@( s1,u1)exp(-theta*(s2-u2)-theta*(s1-u1)).*(s1-s2).^(2*H-2).*(u1-u2).^(2*H-2);

ff5=@(s2,u2)@( s1,u1)exp(-theta*(s2-u2)-theta*(u1-s1)).*(s2-s1).^(2*H-2).*(u2-u1).^(2*H-2);
ff6=@(s2,u2)@( s1,u1)exp(-theta*(s2-u2)-theta*(u1-s1)).*(s2-s1).^(2*H-2).*(u1-u2).^(2*H-2);
ff7=@(s2,u2)@( s1,u1)exp(-theta*(s2-u2)-theta*(u1-s1)).*(s2-s1).^(2*H-2).*(u2-u1).^(2*H-2);
ff8=@(s2,u2)@( s1,u1)exp(-theta*(s2-u2)-theta*(u1-s1)).*(s1-s2).^(2*H-2).*(u2-u1).^(2*H-2);

[email=ff9=@(s2,u2]ff9=@(s2,u2[/email])@( s1,u1)exp(-theta*(u2-s2)-theta*(s1-u1)).*(s2-s1).^(2*H-2).*(u2-u1).^(2*H-2);
ff10=@(s2,u2)@( s1,u1)exp(-theta*(u2-s2)-theta*(s1-u1)).*(s2-s1).^(2*H-2).*(u1-u2).^(2*H-2);
ff11=@(s2,u2)@( s1,u1)exp(-theta*(u2-s2)-theta*(s1-u1)).*(s1-s2).^(2*H-2).*(u2-u1).^(2*H-2);
ff12=@(s2,u2)@( s1,u1)exp(-theta*(u2-s2)-theta*(s1-u1)).*(s1-s2).^(2*H-2).*(u1-u2).^(2*H-2);
ff13=@(s2,u2)@( s1,u1)exp(-theta*(u2-s2)-theta*(u1-s1)).*(s2-s1).^(2*H-2).*(u2-u1).^(2*H-2);
ff14=@(s2,u2)@( s1,u1)exp(-theta*(u2-s2)-theta*(u1-s1)).*(s2-s1).^(2*H-2).*(u1-u2).^(2*H-2);
ff15=@(s2,u2)@( s1,u1)exp(-theta*(u2-s2)-theta*(u1-s1)).*(s2-s1).^(2*H-2).*(u2-u1).^(2*H-2);
ff16=@(s2,u2)@( s1,u1)exp(-theta*(u2-s2)-theta*(u1-s1)).*(s1-s2).^(2*H-2).*(u2-u1).^(2*H-2);

现在请问对应的 Int1……Int16 应该怎么写呢?  谢谢

[ 本帖最后由 ChaChing 于 2010-6-22 23:43 编辑 ]

rocwoods

楼主，你的H为0.78，这样2*H-2为负值，s2-s1可以取到0，这样积分范围内有大量的奇点，积分函数都报错，如果把H改为1.78，2*H-2>0，这样不用把原积分表达式拆成16个，直接按下面按下面方法做即可（T取10会有一些warning，暂时没有细研究，我把T改成1了）

1. 
   
2. theta=1; T=1; H=1.78;
   
3. ff=@(u2,s2)@( u1,s1)exp(-theta*abs(s2-u2)-theta*abs(s1-u1)).*abs(s2-s1).^(2*H-2).*abs(u2-u1).^(2*H-2);
   
4. quad2d(@(u2,s2) arrayfun(@(u2,s2) quad2d(ff(u2,s2),0,T,0,T) ,u2,s2),0,T,0,T)
   
5. 
   
6. ans =
   
7. 
   
8.     0.0248
   

复制代码

weilinhy

谢谢rocwoods的回复
我想请教一个问题：
由于H是赫斯特指数 所以0<H<1。请问这样的话有什么办法么？
非常感谢rocwoods和ChaChing