weilinhy 发表于 2010-4-14 09:58

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


我用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 发表于 2010-4-14 10:48

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

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

weilinhy 发表于 2010-4-14 23:38

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

maigicku 发表于 2010-4-15 09:01

举个例子:
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 发表于 2010-4-15 22:09

谢谢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 发表于 2010-4-16 08:37

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 发表于 2010-4-16 09:49

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 发表于 2010-4-16 11:21

楼主,原被积函数含有绝对值,需要按y和z,x和y的关系拆成四个积分的和,解决代码如下(需要R2009a以上的版本):
clear all
theta=1; T=10; H=0.78;
ff1=@(x) @(y,z) exp(-theta*x-theta*(y-z)).*((x-y)).^(2*H-2).*z.^(2*H-2);
ff2=@(x) @(y,z) exp(-theta*x-theta*(y-z)).*((y-x)).^(2*H-2).*z.^(2*H-2);
ff3=@(x) @(y,z) exp(-theta*x-theta*(z-y)).*((x-y)).^(2*H-2).*z.^(2*H-2);
ff4=@(x) @(y,z) exp(-theta*x-theta*(z-y)).*((y-x)).^(2*H-2).*z.^(2*H-2);
Int1 = quadgk(@(x) arrayfun(@(x) quad2d(ff1(x),0,x,0,@(y)y),x),0,T)
Int2 = quadgk(@(x) arrayfun(@(x) quad2d(ff2(x),x,T,0,@(y)y),x),0,T)
Int3 = quadgk(@(x) arrayfun(@(x) quad2d(ff3(x),0,x,@(y)y,T),x),0,T)
Int4 = quadgk(@(x) arrayfun(@(x) quad2d(ff4(x),x,T,@(y)y,T),x),0,T)
INT = Int1+Int2+Int3+Int4

Int1 =

    1.2617


Int2 =

    4.1067


Int3 =

    1.4130


Int4 =

    3.2933


INT =

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

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

weilinhy 发表于 2010-4-16 11:30

多谢 多谢 rocwoods

maigicku 发表于 2010-4-16 16:32

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 发表于 2010-4-16 17:59

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

weilinhy 发表于 2010-4-17 01:06

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

weilinhy 发表于 2010-6-22 06:55

请教matlab多重积分问题


求解上面积分问题, 参考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);

ff9=@(s2,u2)@( 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 发表于 2010-6-22 10:17

楼主,你的H为0.78,这样2*H-2为负值,s2-s1可以取到0,这样积分范围内有大量的奇点,积分函数都报错,如果把H改为1.78,2*H-2>0,这样不用把原积分表达式拆成16个,直接按下面按下面方法做即可(T取10会有一些warning,暂时没有细研究,我把T改成1了)
theta=1; T=1; H=1.78;
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);
quad2d(@(u2,s2) arrayfun(@(u2,s2) quad2d(ff(u2,s2),0,T,0,T) ,u2,s2),0,T,0,T)

ans =

    0.0248

weilinhy 发表于 2010-6-22 21:19

回复 14楼 rocwoods 的帖子

谢谢rocwoods的回复
我想请教一个问题:
由于H是赫斯特指数 所以0<H<1。请问这样的话有什么办法么?
非常感谢rocwoods和ChaChing
页: [1] 2
查看完整版本: 请教多重积分问题 谢谢大家