声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3460|回复: 17

[编程技巧] 如何用数值方法求以下积分

[复制链接]
发表于 2007-9-1 22:31 | 显示全部楼层 |阅读模式

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

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

x
aa=x*z*sinh(x^2)^0.5 对z在0~0.0001上积分,得关于x的函数。--如果用数值方法?

实在搞不明白了。高手帮忙看看。
clear; clc
syms x z
aa=inline('x.*z.*sinh(x.^2).^0.5','x','z')
quadl(aa,0,0.0001)

提示:
??? Error using ==> inline.feval
Not enough inputs to inline function.

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

Error in ==> h at 14
quadl(aa,0,0.0001)

[ 本帖最后由 ChaChing 于 2010-8-10 15:48 编辑 ]
回复
分享到:

使用道具 举报

发表于 2007-9-2 16:49 | 显示全部楼层
对z积分可以直接用啊, 整个函数都是关于z的一次函数,

aa=x*sinh(x^2)^0.5×5×10^(-9)
 楼主| 发表于 2007-9-2 16:58 | 显示全部楼层
大侠,您是运行的下面这段程序么?
clear
clc
syms x z
aa=inline('x.*z.*sinh(x.^2).^0.5','x','z')
quadl(aa,0,0.0001)
发表于 2007-9-3 08:51 | 显示全部楼层
这个手算就可以了。因为原函数关于z是线性函数。如果非要用数值方法的话这样做就可以了:

  1. aa=@(x) quadl(@(z) x*z*sinh(x^2)^0.5,0,0.0001)
复制代码

评分

1

查看全部评分

 楼主| 发表于 2007-9-3 17:33 | 显示全部楼层

请问,怎样才能把一个表达式代入到一个匿名函数中

谢谢rocwoods大侠.
打算把dpdx的表达式代入到被积函数中,先对z在0到0.0001上积分,再对x在-1到2上积分。运行时总报错。:@(
clear
clc
syms x z
p=(1-x^2)^0.5
dpdx=diff(p)
aa=@(t1,t2)quadl(@(x)quadl(@(z) x.*z.*sinh(dpdx.*z).^0.5,0,0.0001),t1,t2)
bb=aa(-1,2)

[ 本帖最后由 ChaChing 于 2010-8-10 15:52 编辑 ]
发表于 2007-9-3 18:27 | 显示全部楼层

  1. p='(1-x^2)^0.5';
  2. dpdx=diff(p);
  3. fzx=eval([' @(z,x) x.*z.*sinh(',char(vectorize(dpdx)),'.*z).^0.5'])
  4. aa=@(t1,t2) dblquad(fzx,0,0.0001,t1,t2)
复制代码

评分

2

查看全部评分

 楼主| 发表于 2007-9-3 20:53 | 显示全部楼层

回复 #8 rocwoods 的帖子

又学到了不少知识,真的不知道怎么感谢您才好。:handshake
还有一个问题,积分过程中出现了 “被零除” 的提示,是不是说明这个程序有问题啊?

[ 本帖最后由 ChaChing 于 2010-8-11 01:05 编辑 ]
发表于 2007-9-5 21:38 | 显示全部楼层
将积分下限加一个eps试试.
 楼主| 发表于 2007-9-5 21:47 | 显示全部楼层

回复 #11 xjzuo 的帖子

试过了,以前都是NAN+NANi,加eps后,变成NAN了。
发表于 2007-9-5 22:02 | 显示全部楼层
那就是端点奇异,舍弃之,再试试.

另外,我本想针对你的这个问题写一个示例贴的, 总结一下符号求积和数值求积的,但想想其实也很简单,其实以前的一些帖子已经有答案了,只是需要稍微灵活运用一下,所以决定不公布我的算法了.
另:你的这个问题恰好在D之前都有解析表达式,所以符号法和数值法都可用,只是前者需要在D之后,用一次quadl.
根据你贴的公式,结果计算值很小,不知是否在合理范围: ---------(我指的是另一个帖子中的附件word中的公式)
E =1.6162e-006
F =2.5859e-008

[ 本帖最后由 ChaChing 于 2010-8-11 01:07 编辑 ]

评分

1

查看全部评分

 楼主| 发表于 2007-9-5 22:20 | 显示全部楼层
这个问题是我简化了的,赋的值与计算过程中得到的值比较接近,结果比较小是合理的。我写的程序比较繁。应该是端点奇异,被积表达式中的出现了+inf和-inf,但是不知道该怎么舍弃。

刚刚将积分限向里收了一下,没什么效果。

[ 本帖最后由 ChaChing 于 2010-8-11 01:08 编辑 ]
发表于 2007-9-5 23:04 | 显示全部楼层

回复 #5 rocwoods 的帖子

对于这个问题如果用符号积分就非常简单
syms x z;
int(x*z*sinh(x^2)^0.5,z,0,0.0001)

ans =

1/200000000*sinh(x^2)^(1/2)*x
但若要用数值积分,对于有多个参数的数值积分我翻阅了一些参考书,还是不清楚@的用法,烦请高人明示。
发表于 2007-9-5 23:19 | 显示全部楼层


看 xjzuo 版主写的积分、微分问题的几个示例贴(在置顶帖中找)
发表于 2007-9-6 09:27 | 显示全部楼层
注意: 我指的是你另一个帖子中的附件word中的公式.
本贴的例子实在过于简单,而且已有类似的帖子,所以没有回复.
你可以再回去看看你原来的那个帖子中的公式,也就是没有外循环的情形,看看计算结果是否与我的一致.


[ 本帖最后由 xjzuo 于 2007-9-6 09:35 编辑 ]
 楼主| 发表于 2007-9-6 10:02 | 显示全部楼层

回复 #19 xjzuo 的帖子

还是不对呀。
D =

     Inline function:
     D(x) = -1./10000000./(1-x.^2).^(1./2).*x+1000000.*sinh((1000000000000000000.*(1-x.^2).^(1./2).*(exp(1./10000000000000./(1-x.^2).^(1./2).*x)-exp(-1./10000000000000./(1-x.^2).^(1./2).*x))./x+25000000000000000000.*(1-x.^2).^(1./2).*(exp(1./10000000000000./(1-x.^2).^(1./2).*x)+exp(-1./10000000000000./(1-x.^2).^(1./2).*x)-2)./x.*(625000000000000000000000000000000000000.*(1-x.^2).*(exp(1./10000000000000./(1-x.^2).^(1./2).*x)-exp(-1./10000000000000./(1-x.^2).^(1./2).*x)).^2./x.^2-625000000000000000000000000000000000000.*(1-x.^2).*(exp(1./10000000000000./(1-x.^2).^(1./2).*x)+exp(-1./10000000000000./(1-x.^2).^(1./2).*x)-2).^2./x.^2+1./625).^(1./2))./(625000000000000000000000000000000000000.*(1-x.^2).*(exp(1./10000000000000./(1-x.^2).^(1./2).*x)-exp(-1./10000000000000./(1-x.^2).^(1./2).*x)).^2./x.^2-625000000000000000000000000000000000000.*(1-x.^2).*(exp(1./10000000000000./(1-x.^2).^(1./2).*x)+exp(-1./10000000000000./(1-x.^2).^(1./2).*x)-2).^2./x.^2))


E =

      NaN +    NaNi


F =

      NaN +    NaNi

clear ;clc;
syms x y z
f=(1-x.^2).^(1/2);
dfdx=diff(f);
A=int(5e13*cosh(dfdx*z/1e6),z,0,1e-7);
B=int(5e13*sinh(dfdx*z/1e6),z,0,1e-7);
C=(0.04*A-B*(A^2-B^2+0.04^2)^0.5)/(A^2-B^2);
D=dfdx*(1e-7)+1e6*sinh(C);
D=inline(D)
E=quadl(D,(-1e-4),1e-4)
F=0.016*E

我怎么算不出那个结果呢?
C有一点特殊,它是A和B的那样函数,而A和B又是关于x的表达式.

[ 本帖最后由 twomao 于 2007-9-6 10:17 编辑 ]
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-9 04:57 , Processed in 0.071651 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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