matlab求定积分之实例说明

更新时间:2023-05-02 14:45:01 阅读量: 实用文档 文档下载

说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。

一、符号积分

符号积分由函数int来实现。该函数的一般调用格式为:

int(s):没有指定积分变量和积分阶数时,系统按findsym函数指示的默认变量对被积函数或符号表达式s求不定积分;

int(s,v):以v为自变量,对被积函数或符号表达式s求不定积分;

int(s,v,a,b):求定积分运算。a,b分别表示定积分的下限和上限。该函数求被积函数在区间[a,b]上的定积分。a和b可以是两个具体的数,也可以是一个符号表达式,还可以是无穷(inf)。当函数f关于变量x在闭区间[a,b]上可积时,函数返回一个定积分结果。当a,b中有一个是inf时,函数返回一个广义积分。当a,b中有一个符号表达式时,函数返回一个符号函数。

例:

求函数x^2+y^2+z^2的三重积分。内积分上下限都是函数,对z积分下限是sqrt(x*y),积分上限是x^2*y;对y积分下限是sqrt(x),积分上限是x^2;对x的积分下限1,上限是2,求解如下:

>>syms x y z %定义符号变量

>>F2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) %注意定积分的书写格式

F2 =

1610027357/6563700-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2 ^(3/4) %给出有理数解

>>VF2=vpa(F2) %给出默认精度的数值解

VF2 =

224.92153573331143159790710032805

二、数值积分

1.数值积分基本原理

求解定积分的数值方法多种多样,如简单的梯形法、辛普生(Simpson)?法、牛顿-柯特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间[a,b]分成n个子区间[xi,xi+1],i=1,2,…,n,其中x1=a,xn+1=b。这样求定积分问题就分解为求和问题。

2.数值积分的实现方法

基于变步长辛普生法,MATLAB给出了quad函数来求定积分。该函数的调用格式为:

[I,n]=quad('fname',a,b,tol,trace)

基于变步长、牛顿-柯特斯(Newton-Cotes)法,MATLAB给出了quadl函数来求定积分。该函数的调用格式为:

[I,n]=quadl('fname',a,b,tol,trace)

其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度,缺省时取tol=0.001。trace控制是否展现积分过程,若取非0则展现积分过程,取0则不展现,缺省时取trace=0。返回参数I即定积分值,n为被积函数的调用次数。

例:

求函数'exp(-x*x)的定积分,积分下限为0,积分上限为1。

>>fun=inline('exp(-x.*x)','x'); %用内联函数定义被积函数fname

1 / 1

>>Isim=quad(fun,0,1) %辛普生法

Isim =

0.746824180726425

IL=quadl(fun,0,1) %牛顿-柯特斯法

IL =

0.746824133988447

1 / 1

三、梯形法求向量积分

trapz(x,y)—梯形法沿列方向求函数Y关于自变量X的积分(向量形式,数值方法)。

>>d=0.001;

>>x=0:d:1;

>>S=d*trapz(exp(-x.^2))

S=

0.7468

或:

>>format long g

>>x=0:0.001:1; %x向量,也可以是不等间距

>>y=exp(-x.^2); %y向量,也可以不是由已知函数生成的向量

>>S=trapz(x,y); %求向量积分

S =

0.746824071499185

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

int是符号解,无任何误差,唯一问题是计算速度;quad是数值解,有计算精度限制,优点是总是能有一定的速度,即总能在一定时间内给出一个一定精度的解。[FROM: 58.192.116.*]

对于y=exp(-(x.^2+x+1)/(1+x)),被积函数之原函数无"封闭解析表达式",符号计算无法解题,这是符号计算有限性,结果如下:

>> syms x

>>y=exp(-(x.^2+x+1)/(1+x))

>>s=int(y,x,0,inf)

y =

exp((-x^2-x-1)/(1+x))

Warning: Explicit integral could not be found.

>> In sym.int at 58

s =

int(exp((-x^2-x-1)/(1+x)),x = 0 .. Inf)

只有通过数值计算解法

>> dx=0.05; %采样间隔

>>x=0:dx:1000; %数值计算适合于有限区间上,取有限个采样点,只要终值足够大,精度不受影响

>>y=exp(-(x.^2+x+1)./(1+x));

>>S=dx*cumtrapz(y); %计算区间内曲线下图形面积,为小矩形面积累加得

>>S(end)

ans =

0.5641 %所求定积分值

或进行编程,积分上限人工输入,程序如下:

%表达式保存为函数文件

function y=fxy(x)

y=exp(-(x.^2+x+1)./(1+x)); % save fxy.m

% main --------主程序

clear,clc

h=.001;p=0;a=0;

R=input('请输入积分上限,R=')

while a

p=p+(fxy(a)+fxy(a+h))*h/2;

a=a+h;

end

p=vpa(p,10)

运行主程序后得到结果:

请输入积分上限,R=1000

R =

1000

p =

.5641346055

其它结果如下:

0-1: int=.3067601686

0-2: int=.4599633159

0-5: int=.5583068217

0-10: int=.5640928975

0-100: int=.5641346055

0-1000: int=.5641346055

[FROM: 211.65.33.*]

在积分函数中,

sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi*z/9);已知变量e1,e2,e3,n1,n2,n3通过函数参数输入,如果直接用inline或字符串的形式,则表达式中的未知数有9个,分别是e1,e2,e3,n1,n2,n3,x,y,z。而用匿名函数时,已知变量e1,e2,e3,n1,n2,n3就会以常数看待,未知数就只有x,y,z 了,可以求三重积分了。

完整函数程序:

function Fn(n1,n2,n3)

if n1==0

e1=1;

else if n1>0

e1=2;

end

end

if n2==0

e2=1;

else if n2>0

e2=2;

end

end

if n3==0

e3=1;

else if n3>0

e3=2;

end

end

F=@(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi* z/9);

S=triplequad(F,-6,6,-5.5,5.5,-4.5,4.5) %求三重数值积分

将以上代码保存为Fn.m程序文件,即m文件,然后运行:

>> Fn(1,1,1)

S =

866.9655

[FROM: 211.65.33.*]

三重积分请用三重积分函数triplequad,与三个积分上下限对应,即x=triplequad(F,-6,6,-5.5,5.5,-4.5,4.5)

其中被积函数F用"匿名函数"来表达,即

F=@(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi*z/9);

如果直接用inline或字符串的形式,则表达式中的未知数有9个,分别是e1,e2,e3,n1,n2,n3,x,y,z。而用匿名函数时,已知变量e1,e2,e3,n1,n2,n3就会以常数看待,未知数就只有x,y,z了。

完整函数程序:

function Fn(n1,n2,n3)

if n1==0

e1=1;

else if n1>0

e1=2;

end

end

if n2==0

e2=1;

else if n2>0

e2=2;

end

end

if n3==0

e3=1;

else if n3>0

e3=2;

end

end

F=@(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi*z/9);

x=triplequad(F,-6,6,-5.5,5.5,-4.5,4.5)

>> Fn(1,1,1)

x =

866.9655

[FROM: 58.192.116.*]

(注:本资料素材和资料部分来自网络,仅供参考。请预览后才下载,期待您的好评与关注!)

本文来源:https://www.bwwdw.com/article/2rse.html

Top