quad8是matlab中调用那个,Matlab 数值积分
MATLAB数值积分与微分
8.1 数值积分
8.1.1 数值积分基本原理
求解定积分的数值方法多种多样,如简单的梯形法、辛普生(Simpson)法、牛顿-柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间[a,b]分成n个子区间[xi,xi+1],i=1,2,…,n,其中x1=a,xn+1=b。这样求定积分问题就分解为求和问题。
8.1.2 数值积分的实现方法
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函数,从而保证能以更高的效率求出所需的定积分值。
例8-2 求定积分。
(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
例8-3 分别用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)。
例8-4 用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完全相同。
例8-5 计算二重定积分
(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
8.2 数值微分
8.2.1 数值差分与差商
8.2.2 数值微分的实现
在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,按行计算差分。
例8-6 生成以向量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)求数值导数
,diff是后一个数减去前一个数。在最后面加了一个点3.01->才能求出x=3的微分值。n+1个点中间有n段线段
gx=g(x); %求函数f的导函数g在假设点的导数
plot(x,dpx,x,dx,'.',x,gx,'-'); %作图
一、Z = trapz(X,Y,dim)
梯形数值积分,通过已知参数x,y按dim维使用梯形公式进行积分。
例1 计算int(sin(x),0,pi)
>>x=0:pi/100:2*pi;
>>y=sin(x);
>>z=trapz(x,y)%或者说使用z =
pi/100*trapz(y)
z =
1.0300e-017
>>z = pi/100*trapz(y)
二、[q,fcnt]=
quad(fun,a,b,tol,trace,p1,p2...)
自适应simpson公式数值积分,适用于精度要求低,被积函数平滑性较差的数值积分。
注意事项:
1.被积函数fun必须是函数句柄;
2.积分限[a,b]必须是有限的,因此不能为inf;
3.p1为其他需要传递的参数,一般是数值。
可能警告:
1.'Minimum step size reached'
意味着子区间的长度与计算机舍入误差相当,无法继续计算了。原因可能是有不可积的奇点;
2.'Maximum function count exceeded'
意味着积分递归计算超过了10000次。原因可能是有不可积的奇点;
3.'Infinite or Not-a-Number function value encountered'
意味着在积分计算时,区间内出现了浮点数溢出或者被零除。
例2 计算积分1/(x^3-2*x-p),其中参数p=5,积分区间为[0,2]。
>>F = @(x,n)1./(x.^3-2*x-n);
>>Q = quad(@(x)F(x,5),0,2)%或者使用
quad(F,0,2,[],[],5)效果是一样的,只是前者使用的函数嵌套。
Q =
-0.4605
>>quad(F,0,2,[],[],5)
ans =
-0.4605
三、[q,fcnt] =
quadl(fun,a,b,tol,trace,p1,p2...)
自适应Lobatto数值积分,适用于精度要求高,被积函数曲线比较平滑的数值积分。
注意事项:
同quad
可能警告:
同quad
例3 计算积分1/(x^3-2*x-p),其中参数p=5,积分区间为[0,2]。
>>F=@(x,p)1./(x.^3-2*x-p);
>>Q = quadl(F,0,2,[],[],5)%或者Q =
quadl(@(x)F(x,5),0,2)。
Q =
-0.4605
四、[q,errbnd] =
quadgk(fun,a,b,param1,val1,param2,val2,...)
自适应Gauss-Kronrod数值积分,适用于高精度和震荡数值积分,支持无穷区间,并且能够处理端点包含奇点的情况,同时还支持沿着不连续函数积分,复数域线性路径的围道积分法。
注意事项:
1.积分限[a,b]可以是[-inf,inf],但必须快速衰减;
2.被积函数在端点可以有奇点,如果区间内部有奇点,将以奇点区间划分成多个,也就是说奇点只能出现在端点上;
3.被积函数可以剧烈震荡;
4.可以计算不连续积分,此时需要用到'Waypoints'参数,'Waypoints'中的点必须严格单调;
5.可以计算围道积分,此时需要用到'Waypoints'参数,并且为复数,各点之间使用直线连接;
6.param,val为函数的其它控制参数,比如上面的'waypoints'就是,具体看帮助。
出现错误:
1.'Reached the limit on the maximum number of intervals in
use'
2.'Infinite or Not-a-Number function value encountered'
例4 计算有奇点积分int(exp(x)*log(x),0,1)。
>>F=@(x)exp(x).*log(x);%奇点必须在端点上,否则请先进行区间划分。
>>Q = quadgk(F,0,1)
Q =
-1.3179
例5 计算半无限震荡积分int(x^5*exp(-x)*sin(x),0,inf)。
>>F=@(x)x.^5.*exp(-x).*sin(x);
>>fplot(F,[0,100])%绘图,看看函数的图形
>>[q,errbnd] =
quadgk(F,0,inf,'RelTol',1e-8,'AbsTol',1e-12)%积分限中可以有inf,但必须快速收敛
q =
-15.0000
errbnd =
9.4386e-009
例6
计算不连续积分,积分函数为f(x)=x^5*exp(-x)*sin(x),但是人为定义f(2)=1000,f(5)=-100,积分区间为[1
10]。
>>F=@(x)x.^5.*exp(-x).*sin(x);
>>[q,errbnd] =
quadgk(F,1,10,'Waypoints',[2 5])%显然2,5为间断点
q =
-10.9408
errbnd =
3.2296e-014
例7 计算围道积分,在复数域内,积分函数1/(2*z-1),积分路径为由[-1-i 1-i 1+i -1+i
-1-i]围成的矩形边框。
>>Waypoints=[-1-i 1-i 1+i -1+i
-1-i];
>>plot(Waypoints);%绘制积分路径
>>xlabel('Real axis');ylabel('Image
axis');axis([-1.5 1.5 -1.5 1.5]);grid on;
>>Q = quadgk(@(z)1./(2*z -
1),-1-i,-1-i,'Waypoints',[1-i,1+i,-1+i])%注意各点间使用直线连接
ans =
0.0000 + 3.1416i
>> quadgk(@(z)1./(2*z -
1),-1-i,-1-i,'Waypoints',Waypoints)%使用这个的效果也是一样的,就是说始末点可以随便包不包含在Waypoints中。
ans =
0.0000 + 3.1416i
五、[Q,fcnt] = quadv(fun,a,b,tol,trace)
矢量化自适应simpson数值积分。
注意事项:
1.该函将quad函数矢量化了,就是一次可以计算多个积分;
2.所有的要求完全与quad相同。
例8 计算下面积分,分别计算n=1,2...,5时的5个积分值,被积函数1/(n+x),积分限为[0,1]。
>>for k =
1:5, Qs(k) = quadv(@(x)1/(k+x),0,1);end;Qs
Qs =
0.6931 0.4055 0.2877 0.2231 0.1823
>>F=@(x,n)1./((1:n)+x);%定义被积函数
>>quadv(@(x)F(x,5),0,1)%我们可以完全使用quadv函数替换上面循环语句的,建议使用后者
ans =
0.6931 0.4055 0.2877 0.2231 0.1823
六、q =
dblquad(fun,xmin,xmax,ymin,ymax,tol,method)
矩形区域二重数值积分,一般区域二重积分参见NIT(数值积分工具箱)的quad2dggen函数。
例9 计算下面二重积分
>>F = @(x,y)y*sin(x)+x*cos(y);
>Q = dblquad(F,pi,2*pi,0,pi)
Q =
-9.8696
七、q=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)
长方体区域三重数值积分,注意此时没有一般区域的三重积分。
例10 计算下面三重积分。
>>F =
@(x,y,z)y*sin(x)+z*cos(x);
>>Q = triplequad(F,0,pi,0,1,-1,1)
Q =
2.0000
八、超维长方体区域多重积分
quadndg:NIT工具箱函数,可以解决多重超维长方体边界的定积分问题,但没有现成的一般积分区域求解函数。
下面总结下:
(1)quad:采用自适应变步长simpson方法,速度和精度都是最差的,建议不要使用;
(2)quad8:使用8阶Newton-Cotes算法,精度和速度均优于quad,但在目前版本下已被取消;
(3)quadl:采用lobbato算法,精度和速度均较好,建议全部使用该函数;
(4)quadg:NIT(数值积分)工具箱函数,效率最高,但该工具箱需要另外下载;
(5)quadv:quad的矢量化函数,可以同时计算多个积分;
(6)quadgk:很有用的函数,功能在Matlab中最强大;
(7)quad2dggen:一般区域二重积分,效率很好,需要NIT支持;
(8)dblquad:长方形区域二重积分;
(9)triplequadL:长方体区域三重积分;
(10)quadndg:超维长方体区域积分,需要NIT支持。
int的积分可以是定积分,也可以是不定积分(即有没有积分上下限都可以积),可以得到解析的解,比如你对x^2积分,得到的结果是1/3*x^3,这是通过解析的方法来解的。如果int(x^2,x,1,2)得到的结果是7/3
;
quad是数值积分,它只能是定积分(就是有积分上下限的积分),它是通过simpson数值积分来求得的(并不是通过解析的方法得到解析解,再将上下限代入,而是用小梯形的面积求和得到的)。如果f=inline(‘x.^2′);quad(f,1,2)得到的结果是2.333333,这个数并不是7/3
;
最新心得:
看一本书上介绍quad积分时,是创建了一个子函数形式,如
function f=quadl(x)
f=x.^2;
Q=quad(‘quadl’,0,2)
结果Q =
2.6667
如果函数中有一个已知变量如a的话,如
function f=quadl(x)
a=2;
f=a+x.^2;
Q=quad(‘quadl’,0,2)
结果Q =
6.6667
当用使用inline函数的时候可以避免调用子函数的麻烦,直接把这个功能集成于总程序,如
f=inline(‘x.^2′);quad(f,1,2)
但是当函数为
a=2;
f=inline(‘a+x.^2′);
quad(f,1,2)
计算就会出错,说明inline中不能带已知的字母。
但是很多时候,变量a是循环变化的,这样就导致这种调用子函数的方法非常不好用,
a不能及时改变,下面的方法可以解决这个问题:
a=2;
f=@(x)(a+x.^2);
Q=quad(f,0,2)
结果Q =
6.6667 正确
用@来表达函数要比inline应用更广,在循环中应用更有利!
quad8是matlab中调用那个,Matlab 数值积分相关推荐
- 怎么用matlab编写quad8算法,quad8 是 matlab 中调用( )公式的命令
quad8 是 matlab 中调用( )公式的命令 更多相关问题 [判断题]对于糖尿病患者来说,心脏病.卒中和猝死的风险减少一半. [单选题]孔子曰:人皆可为尧舜.只要我们从小形成良好的生活态度.人 ...
- matlab中调用java代码_Matlab中调用第三方Java代码
在Java中采用Matlab JA Builder可以实现调用m文件,采用这样的方式,可在Matlab的M文件中,直接调用Java类.这种方式可以表示为Java--> Matlab( m, Ja ...
- 在Matlab中调用Mathematica的函数
Matlab与Mathematica都是目前比较好的数学软件,个人感觉Matlab比较偏应用,工程计算方面,而Mathematica比较适合数学推理,,理论比较重,个人最喜欢就是它的二维输入格式,比较 ...
- matlab中调用cplex 以及使用 Yalmip 工具箱
matlab中调用cplex 以及使用 Yalmip 工具箱_天天向上的专栏-CSDN博客_matlab调用cplex CPLEX 在matlab中实现的一个例子-RCPSP_starry0001的博 ...
- m 文件 dll matlab 中调用_如何在matlab中调用python程序
现在python很火,很多代码都是python写的,如果你和我一样,习惯了使用matlab,还想在matlab中调用Python的代码,应该怎么办呢?其中一条思路:首先在matlab中调用系统脚本命令 ...
- php调用python绘图程序_如何在matlab中调用python程序
现在python很火,很多代码都是python写的,如果你和我一样,习惯了使用matlab,还想在matlab中调用Python的代码,应该怎么办呢?其中一条思路:首先在matlab中调用系统脚本命令 ...
- 将keras训练好的神经网络模型生成独立于python环境运行的.exe文件并在matlab中调用(上篇)
keras是深度学习简单易上手的python第三方库.在实际应用中,常常需要将keras训练好的深度学习模型生成独立于python环境的文件,供不会python的人或未安装python环境的设备调用. ...
- matlab实时摄像头处理,在MATLAB中调用摄像头实时地显示图像
<在MATLAB中调用摄像头实时地显示图像>由会员分享,可在线阅读,更多相关<在MATLAB中调用摄像头实时地显示图像(1页珍藏版)>请在人人文库网上搜索. 1.clc;clf ...
- 如何在MATLAB中调用(运行)“用Python写成的函数或脚本”
如何在MATLAB中调用"用Python写成的函数或脚本",首先要确保MATLAB知道咱们的Python解释器的位置在哪里. 如果安装了Python的时候把Python加入了系统环 ...
最新文章
- 关系型数据库之Mysql备份(五)
- Windows server 2012 搭建×××图文教程(一)安装×××相关服务
- 5G NR — 关键技术
- 网站优化之如何筛选更正确高质量的关键词?
- Spring学习8之静态代理
- 1.5 编程基础之循环控制 44 第n小的质数
- GitHub上如何创建文件夹
- (十)OpenCV相机标定
- 我们自研的那些Devops工具
- 摩托罗拉里程碑2刷android 2.3.4
- 这几个小众软件能支持英文文档翻译成中文
- Python 实现文字聊天室-功能拓展
- 灵雀云ACP成功通过金融信创生态实验室适配验证
- linux系统vi替换命令,linux:vi 替换命令
- linux的标准输入输出
- 源代码静态检测分析技术浅析
- 发生系统错误 5。拒绝访问。
- 时间序列预测分析(2)How to Develop a Skillful Machine Learning Time Series Forecasting Model
- web前端工程师简历
- Airship再下一城,因为戴尔和AtT的合作
热门文章
- python中的any与all函数的区别
- Python中九种格式化输出方法,你都知道吗?
- 10行Python代码自动清理电脑内重复文件
- html 显示搜索结果,搜索结果高亮显示(不改变html标签)
- 【植物大战僵尸2】算法 笔记
- 如何在C++中调用python程序?
- Intel Realsense D435 在C/C++中表示的frame_set就是python中的frames?【wait_for_frames()】
- Intel Realsense D435 pipeline对象指什么?
- PHP中substr截取中文乱码解决方案
- python的读取纯文本文件的几种模式