声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4284|回复: 22

[综合讨论] 请教多重积分问题 谢谢大家

[复制链接]
发表于 2010-4-14 09:58 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
积分.jpg
我用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)如果不行,数值积分求解应该怎么求呢!
谢谢大家

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2010-4-14 10:48 | 显示全部楼层
关于出错信息的解释与改正请参考下帖:
常见的程序出错问题整理

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

评分

1

查看全部评分

 楼主| 发表于 2010-4-14 23:38 | 显示全部楼层
匿名函数怎么弄呢?拜谢阿!刚接触matlab  很多不懂  谢谢了
发表于 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的书。

评分

1

查看全部评分

 楼主| 发表于 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.
帮忙看看 拜谢了
发表于 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 编辑 ]

评分

1

查看全部评分

 楼主| 发表于 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 编辑 ]
发表于 2010-4-16 11:21 | 显示全部楼层
楼主,原被积函数含有绝对值,需要按y和z,x和y的关系拆成四个积分的和,解决代码如下(需要R2009a以上的版本):

  1. clear all
  2. theta=1; T=10; H=0.78;
  3. ff1=@(x) @(y,z) exp(-theta*x-theta*(y-z)).*((x-y)).^(2*H-2).*z.^(2*H-2);
  4. ff2=@(x) @(y,z) exp(-theta*x-theta*(y-z)).*((y-x)).^(2*H-2).*z.^(2*H-2);
  5. ff3=@(x) @(y,z) exp(-theta*x-theta*(z-y)).*((x-y)).^(2*H-2).*z.^(2*H-2);
  6. ff4=@(x) @(y,z) exp(-theta*x-theta*(z-y)).*((y-x)).^(2*H-2).*z.^(2*H-2);
  7. Int1 = quadgk(@(x) arrayfun(@(x) quad2d(ff1(x),0,x,0,@(y)y),x),0,T)
  8. Int2 = quadgk(@(x) arrayfun(@(x) quad2d(ff2(x),x,T,0,@(y)y),x),0,T)
  9. Int3 = quadgk(@(x) arrayfun(@(x) quad2d(ff3(x),0,x,@(y)y,T),x),0,T)
  10. Int4 = quadgk(@(x) arrayfun(@(x) quad2d(ff4(x),x,T,@(y)y,T),x),0,T)
  11. INT = Int1+Int2+Int3+Int4

  12. Int1 =

  13.     1.2617


  14. Int2 =

  15.     4.1067


  16. Int3 =

  17.     1.4130


  18. Int4 =

  19.     3.2933


  20. INT =

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

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

评分

1

查看全部评分

 楼主| 发表于 2010-4-16 11:30 | 显示全部楼层
多谢 多谢 rocwoods
发表于 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
发表于 2010-4-16 17:59 | 显示全部楼层
呵呵,别叫大神,大家都是MATLAB学习者:)
那种用法不是R2009a特有,这是双重匿名函数的用法。7.0以上的版本都支持。主要是quad2d函数是2009a才有的,上面用到这个函数了。不用这个函数也可以,就是稍微麻烦些,要用两次arrayfun,书里面有介绍。
拟合的问题可以问下xiezhh,谢老师很在行的。呵呵
 楼主| 发表于 2010-4-17 01:06 | 显示全部楼层
给您添麻烦了,谢谢您的回答
 楼主| 发表于 2010-6-22 06:55 | 显示全部楼层

请教matlab多重积分问题

积分4.jpg
求解上面积分问题, 参考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 编辑 ]
发表于 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了)

  1. theta=1; T=1; H=1.78;
  2. 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);
  3. quad2d(@(u2,s2) arrayfun(@(u2,s2) quad2d(ff(u2,s2),0,T,0,T) ,u2,s2),0,T,0,T)

  4. ans =

  5.     0.0248
复制代码

评分

1

查看全部评分

 楼主| 发表于 2010-6-22 21:19 | 显示全部楼层

回复 14楼 rocwoods 的帖子

谢谢rocwoods的回复
我想请教一个问题:
由于H是赫斯特指数 所以0<H<1。请问这样的话有什么办法么?
非常感谢rocwoods和ChaChing
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-24 16:37 , Processed in 0.080768 second(s), 31 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表