如何用数值方法求以下积分

twomao

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 编辑 ]

pengweicai

对z积分可以直接用啊， 整个函数都是关于z的一次函数，

aa=x*sinh(x^2)^0.5×5×10^(-9)

twomao

大侠，您是运行的下面这段程序么？
clear
clc
syms x z
aa=inline('x.*z.*sinh(x.^2).^0.5','x','z')
quadl(aa,0,0.0001)

rocwoods

这个手算就可以了。因为原函数关于z是线性函数。如果非要用数值方法的话这样做就可以了：

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

复制代码

twomao

谢谢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 编辑 ]

rocwoods

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

复制代码

twomao

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

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

xjzuo

将积分下限加一个eps试试．

twomao

试过了，以前都是NAN+NANi，加eps后，变成NAN了。

xjzuo

那就是端点奇异,舍弃之,再试试.

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

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

twomao

这个问题是我简化了的，赋的值与计算过程中得到的值比较接近，结果比较小是合理的。我写的程序比较繁。应该是端点奇异，被积表达式中的出现了+inf和-inf，但是不知道该怎么舍弃。

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

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

xyzhou1234

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

ans =

1/200000000*sinh(x^2)^(1/2)*x
但若要用数值积分，对于有多个参数的数值积分我翻阅了一些参考书，还是不清楚@的用法，烦请高人明示。

eight

原帖由 xyzhou1234 于 2007-9-5 23:04 发表

                               
登录/注册后可看大图

但若要用数值积分，对于有多个参数的数值积分我翻阅了一些参考书，还是不清楚@的用法，烦请高人明示

看 xjzuo 版主写的积分、微分问题的几个示例贴（在置顶帖中找）

xjzuo

注意: 我指的是你另一个帖子中的附件word中的公式.
本贴的例子实在过于简单,而且已有类似的帖子,所以没有回复.
你可以再回去看看你原来的那个帖子中的公式,也就是没有外循环的情形,看看计算结果是否与我的一致.


[ 本帖最后由 xjzuo 于 2007-9-6 09:35 编辑 ]

twomao

还是不对呀。
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 编辑 ]