一.数值积分的实现方法

  1．变步长辛普生法
  基于变步长辛普生法，MATLAB给出了quad函数来求定积分。该函数的调用格式为：
  [I,n]=quad('fname',a,b,tol,trace)
  其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。

  例8-1 求定积分。
  (1) 建立被积函数文件fesin.m。
  function f=fesin(x)
  f=exp(-0.5*x).*sin(x+pi/6);
  (2) 调用数值积分函数quad求定积分。
  [S,n]=quad('fesin',0,3*pi)
  S = 0.9008
  n = 77

  2．牛顿－柯特斯法
  基于牛顿－柯特斯法，MATLAB给出了quad8函数来求定积分。该函数的调用格式为：
  [I,n]=quad8('fname',a,b,tol,trace)
  其中参数的含义和quad函数相似，只是tol的缺省值取10-6。该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。
  (1) 被积函数文件fx.m。
  function f=fx(x)
  f=x.*sin(x)./(1+cos(x).*cos(x));
  (2) 调用函数quad8求定积分。
  I=quad8('fx',0,pi)
  I = 2.4674
  分别用quad函数和quad8函数求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。
  调用函数quad求定积分：
  format long;
  fx=inline('exp(-x)');
  [I,n]=quad(fx,1,2.5,1e-10)
  I = 0.28579444254766
  n = 65
  调用函数quad8求定积分：
  format long;
  fx=inline('exp(-x)');
  [I,n]=quad8(fx,1,2.5,1e-10)
  I = 0.28579444254754
  n = 33

  3．被积函数由一个表格定义
  在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X,Y定义函数关系Y=f(X)。
  用trapz函数计算定积分。
  命令如下：
  x=1:0.01:2.5;
        Y=exp(-X); %生成函数关系数据向量
  trapz(X,Y)
  ans = 0.28579682416393

  8.1.3 二重定积分的数值求解
  使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：
  I=dblquad(f,a,b,c,d,tol,trace)
  该函数求f(x,y)在[a,b]×[c,d]区域上的二重定积分。参数tol，trace的用法与函数quad完全相
同。
  计算二重定积分
  (1) 建立一个函数文件fxy.m：
  function f=fxy(x,y)
  global ki;
  ki=ki+1; %ki用于统计被积函数的调用次数
  f=exp(-x.^2/2).*sin(x.^2+y);
  (2) 调用dblquad函数求解。
  global ki;ki=0;
  I=dblquad('fxy',-2,2,-1,1)
  ki
  I = 1.57449318974494
  ki = 1038

  二.数值微分
　
  数值微分的实现
  在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为：
  DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,…,n-1。
  DX=diff(X,n)：计算X的n阶向前差分。例如，diff(X,2)=diff(diff(X))。
  DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。
  生成以向量V=[1,2,3,4,5,6]为基础的范得蒙矩阵，按列进行差分运算。
  命令如下：
  V=vander(1:6)
  DV=diff(V) %计算V的一阶差分

  例8-7 用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f'(x)的图像。
  程序如下：
  f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');
  g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5');
  x=-3:0.01:3;
  p=polyfit(x,f(x),5); %用5次多项式p拟合f(x)
  dp=polyder(p); %对拟合多项式p求导数dp
  dpx=polyval(dp,x); %求dp在假设点的函数值
  dx=diff(f([x,3.01]))/0.01; %直接对f(x)求数值导数
  gx=g(x); %求函数f的导函数g在假设点的导数
  plot(x,dpx,x,dx,'.',x,gx,'-'); %作图
Matlab数值积分的一些做法
Filed under: 未分类 — franz @ 9:31 pm
今天是元宵节，突然来讲Matlab的确很奇怪，连我自己也这样感觉．
由于Matlab对于计算过程的细节要求定义详细，往往让人觉得使用不方便．与同样很流行的Mathematic相比，Matlab更象是程序，而不是象Mathematic那么直观的写作业．
Matlab同样提供符号积分（不定积分）和数值积分（定积分）两大功能．符号积分使用int命令，配合之前的函数定义使用．比如：
a=…;
b=…;
function f = fun(x)
f=a*x^2+b;
将此文件寸为一个单独m文件，再在主程序中调用即可：
F = int(fun,c,d);
c和d为积分上下限，通常为符号，可使用syms c;语句定义。在完成符号积分之后可以对其附值，则就完成了数值积分的任务．
有一点需要注意的是，通过这样过程求的数值，其格式为符号格式’sym’，不能做图，不能和数值进行运算．处理方法是：
vpa (F);　求得３２位符号近似解，再用double ( vpa (F) );将其转换成Matlab默认的双精度位数就可以．注意，直接做 double (F) 行不通，Matlab会提示你＂Conversion from sym to double is not possible＂，不知道＂not possbile"和"impossible"有什么区别，呵呵．
符号积分中一个最大的问题在于存在不可积的函数，比如十分常用的Boltzman分布，或者叫Arrhenius公式：exp(- q*V/ (k *T ))．
插句话，我一直以为Arrhenius是德国人，知道前几天在和老师讨论中讲起，我老师很自豪的对我说，有个著名的瑞典化学家，得过Noble Prize，叫Arrhenius你知道不知道，才发现这个小问题，３滴汗，搞的我老师也很有挫败感．．． diff函数
对于不可积的函数，使用 int 命令之后，得到的表达式中会含有　Ei　项，其本身也是一个函数，以你给定的符号积分上下限为变量．在含有　Ei　项后，对符号上下限附值亦无法得到数值积分解，因为Ei　项也需要解．解Ei　项的方法如下：
result　=　str2num　(maple(‘evalf(Ei(a,b))’)); 
％此语句可以直接使用，其中的a和b就是你所要解数值解的积分上下限，并且只能是具体
数值，不能为符号．
若是积分上限或者下限是一个变化的值　（比如今天我做的一个很简单的计算就是这样的情况），则可以使用以下的方法：
>> result=maple(‘evalf’,'(Ei(1,y))’)
result =Ei(1,y)
>> y=2
y = 2
>> result=subs(result)
result =Ei(1,2)
>> result=vpa(result)
result =.48900510708061119567239835228050e-1
>> result=str2num(maple(‘evalf’,'(Ei(1,2))’))
result = 0.0489
比如其中的y就可以是一个变化的值，写成一个循环，从而计算相应的积分值．
 
Matlab也直接提供两种主要的数值积分函数：quad　和　quad8．quad是变步长辛普生法，quad8是牛顿－柯特斯法，以我今天的例子持续来看，相差很小．
对被积分函数的定义同上，另外，还有一种办法，可以不用另外写一个m文件再调用，省下些小麻烦．可使用 inline 语句：
fxy = inline ( ‘exp(-a*x*y / b)’, ‘x’,'y’);
Matlab将字符串中的表达式录入内存，准备之后使用．这个语句有一个很大的缺点是，表达式中，除了变量之外，其他数值（如上式中的a和b)不能使用字母等符号表示，而必须为数值，即使你已经在之前定义过，也不行，因为　inline语句将引号中的表达式录为符号格
式，其中任何的字母或者字母组合，都会被认为是变量，即使你在语句之后只写了真正的变量，程序还是会全部误读．这就在做积分时候带来错误，明明是常量的a和b，Matlab还是要求你给他们一个积分空间进行积分，从而出错．

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。