如何用数值方法求以下积分
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 编辑 ] 对z积分可以直接用啊， 整个函数都是关于z的一次函数，
aa=x*sinh(x^2)^0.5×5×10^(-9) 大侠，您是运行的下面这段程序么？
clear
clc
syms x z
aa=inline('x.*z.*sinh(x.^2).^0.5','x','z')
quadl(aa,0,0.0001) 这个手算就可以了。因为原函数关于z是线性函数。如果非要用数值方法的话这样做就可以了：
aa=@(x) quadl(@(z) x*z*sinh(x^2)^0.5,0,0.0001)
请问,怎样才能把一个表达式代入到一个匿名函数中
谢谢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 编辑 ]
p='(1-x^2)^0.5';
dpdx=diff(p);
fzx=eval([' @(z,x) x.*z.*sinh(',char(vectorize(dpdx)),'.*z).^0.5'])
aa=@(t1,t2) dblquad(fzx,0,0.0001,t1,t2)
回复 #8 rocwoods 的帖子
又学到了不少知识，真的不知道怎么感谢您才好。:handshake还有一个问题，积分过程中出现了 “被零除” 的提示，是不是说明这个程序有问题啊？
[ 本帖最后由 ChaChing 于 2010-8-11 01:05 编辑 ] 将积分下限加一个eps试试．
回复 #11 xjzuo 的帖子
试过了，以前都是NAN+NANi，加eps后，变成NAN了。 那就是端点奇异,舍弃之,再试试.另外,我本想针对你的这个问题写一个示例贴的, 总结一下符号求积和数值求积的，但想想其实也很简单，其实以前的一些帖子已经有答案了，只是需要稍微灵活运用一下，所以决定不公布我的算法了．
另：你的这个问题恰好在D之前都有解析表达式，所以符号法和数值法都可用，只是前者需要在D之后，用一次quadl.
根据你贴的公式，结果计算值很小，不知是否在合理范围： ---------(我指的是另一个帖子中的附件word中的公式)
E =1.6162e-006
F =2.5859e-008
[ 本帖最后由 ChaChing 于 2010-8-11 01:07 编辑 ] 这个问题是我简化了的，赋的值与计算过程中得到的值比较接近，结果比较小是合理的。我写的程序比较繁。应该是端点奇异，被积表达式中的出现了+inf和-inf，但是不知道该怎么舍弃。
刚刚将积分限向里收了一下，没什么效果。
[ 本帖最后由 ChaChing 于 2010-8-11 01:08 编辑 ]
回复 #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
但若要用数值积分，对于有多个参数的数值积分我翻阅了一些参考书，还是不清楚@的用法，烦请高人明示。 原帖由 xyzhou1234 于 2007-9-5 23:04 发表 http://www.chinavib.com/forum/images/common/back.gif
但若要用数值积分，对于有多个参数的数值积分我翻阅了一些参考书，还是不清楚@的用法，烦请高人明示
看 xjzuo 版主写的积分、微分问题的几个示例贴（在置顶帖中找） 注意: 我指的是你另一个帖子中的附件word中的公式.
本贴的例子实在过于简单,而且已有类似的帖子,所以没有回复.
你可以再回去看看你原来的那个帖子中的公式,也就是没有外循环的情形,看看计算结果是否与我的一致.
[ 本帖最后由 xjzuo 于 2007-9-6 09:35 编辑 ]
回复 #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
我怎么算不出那个结果呢？
Ｃ有一点特殊，它是Ａ和Ｂ的那样函数，而Ａ和Ｂ又是关于x的表达式．
[ 本帖最后由 twomao 于 2007-9-6 10:17 编辑 ]
页:
[1]
2