三个参数matlab程序,用matlab求定积分的三个实例代码
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
三个参数matlab程序,⽤matlab求定积分的三个实例代码
⼀、符号积分
符号积分由函数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
>>Isim=quad(fun,0,1) %⾟普⽣法
Isim =
0.746824180726425
IL=quadl(fun,0,1) %⽜顿-柯特斯法
IL =
0.746824133988447
三、梯形法求向量积分
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.*]。