MATLAB计算积分

函数的积分和椭圆的周长1．正弦函数的积分[问题]求正弦函数从0到π的积分y = sin x当x = 0时，积分为0，画出积分的函数曲线。

[数学模型]定积分的结果为ππ00sin d cos 2S x x x ==-=⎰ 不定积分的结果为sin d cos I x x x C ==-+⎰其中C 是积分常量，由初始条件决定。

当x = 0时，积分为I = 0，必有C = 1。

结果为I = -cos x + 1[算法]根据积分的基本概念，将积分区域分为多份，用矩形法求曲线下的近似面积表示积分的近似值1()ni i S f x x ==∆∑矩形法的函数是sum(f)。

用梯形法求曲线下的近似面积表示积分的近似值1101[()()]2n i i i S f x f x x -+==+∆∑梯形法的函数是trapz(f)。

用数值积分的函数是quad 和quadl ，常用使用格式是S = quad(f,a,b)其中，f 表示被积函数，a 表示积分的下限，b 表示积分的下限。

用符号的函数是int ，常用使用格式是S = int(f,a,b)[程序]zqy4_1.m 如下。

%正弦函数的积分clear %清除变量x=linspace(0,pi); %自变量向量dx=x(2); %间隔y=sin(x); %被积函数s1=sum(y)*dx %矩形法积分s2=trapz(y)*dx %梯形法积分f=inline('sin(x)'); %被积的内线函数s3=quad(f,0,pi) %数值定积分s4=int('sin(x)',0,pi) %符号积分sc1=cumsum(y)*dx; %矩形法累积积分(精度稍差)sc2=cumtrapz(y)*dx; %梯形法累积积分figure %创建图形窗口plot(x,-cos(x)+1,x,sc1,'.',x,sc2,'o') %画解析式和矩阵法以及梯形法积分曲线 s=int('sin(x)') %符号积分sc3=subs(s,'x',x); %替换数值求符号积分的值C=-sc3(1) %求积分常数hold on %保持图像plot(x,sc3+C,'c*') %画符号法积分曲线grid on %加网格fs=16; %字体大小xlabel('\itx','FontSize',fs) %横坐标ylabel('\intsin\itx\rmd\itx','FontSize',fs)%纵坐标title('正弦函数的积分','FontSize',fs) %标题legend('解析解','矩形法','梯形法','符号法')%图例zqy4.1图 zqy4.2图2.三角函数和指数的积分[问题]求如下函数的积分y = e ax sin bx其中a = 0.5，b = 2。

积分下限为0。

画出积分的函数曲线。

[数学模型]设 11e sin d sin de {e sin e cos d }ax ax ax ax I bx x bx bx b bx x a a===-⎰⎰⎰ 11{e sin cos de }{e sin [e cos e sin d ]}ax ax ax ax ax b b bx bx bx bx b bx x a a a a =-=-+⎰⎰ 因此不定积分为221e (sin cos )ax I a bx b bx C a b=-++ 当x = 0时，I 应该为零，所以 22b C a b=+因此，从0开始的积分为221e (sin cos )ax I a bx b bx b a b=-++ [算法]可用积分的解析式直接画图，也可利用被积函数通过梯形求和指令cumtrapz 求积分，数值积分指令quad 求积分，还可通过符号积分int 求解。

[程序]zqy4_2cumtrapz.m 如下。

%数值积分和符号积分方法clear %清除变量a=0.5; %指数的常数b=2; %正弦函数的常数dx=0.1; %间隔xm=6; %上限x=0:dx:xm; %自变量向量s1=(exp(a*x).*(-b*cos(b*x)+a*sin(b*x)))/(a^2+b^2);%积分的解析解C=-s1(1); %求积分常数y=exp(a*x).*sin(b*x); %被积函数s2=cumtrapz(y)*dx; %梯形法积分%s2=cumsum(y)*dx; %矩形法求积分figure %创建图形窗口plot(x,s1+C,x,s2,'.') %画积分曲线s=['exp(',num2str(a),'*x).*sin(',num2str(b),'*x)'];%被积分函数字符串f=inline(s); %化为内线函数s3=0; %第1个值for i=2:length(x) %按自变量循环s3=[s3,quad(f,0,x(i))]; %连接积分end %结束循环s=int('exp(a*x)*sin(b*x)') %对x 进行符号积分ss=subs(s,{'a','b'},{a,b}); %替换常数s4=subs(ss,'x',x); %替换向量hold on %保持图像plot(x,s3,'or',x,s4-s4(1),'ch','MarkerSize',16)%画数值积分和符号积分曲线 title(['\ity\rm=e^{',num2str(a),'}\it^x\rmsin',num2str(b),'\itx\rm 的积分'],'FontSize',16)%标题legend('公式法','梯形法','数值法','符号法',2)%加图例3.椭圆周长的计算[问题]椭圆的长半轴为a ，偏心率为e ，求其周长。

[数学模型]椭圆的方程为x = a sin φ，y = b cos φ其微分为d x = a cos φd φ，y = -b sin φd φ弧元为d s =ϕϕ==ϕ=ϕ= 周长为π/4C a ϕ=⎰其中e 是偏心率。

对于一个偏心率，可直接计算上述积分。

第二类完全椭圆积分定义为π0E()k x =⎰周长可表示为 4E()C a e =[算法]取长半轴为a 为周长的单位，周长可表示为π/*04C C a ϕ==⎰ 可利用数值积分指令quadl 直接对上式积分，也可利用符号积分指令int 对上式积分。

将参数e 2向量和积分变量φ化为矩阵，则能利用梯形求和指令trapz 求周长。

用MATLAB 的函数ellipke 可直接计算第二类完全椭圆积分之值，不过，第二类完全椭圆积分定义为π0E()m x =⎰因此椭圆周长为 24E()C a e =[程序]zqy4_3_1.m 如下。

%椭圆的周长clear %清除变量e=0:0.1:1; %偏心率c1=[]; %周长向量置空for i=1:length(e) %按偏心率循环s=['sqrt(1-',num2str(e(i)^2),'*sin(x).^2)'];%被积函数字符串f=inline(s); %被积内线函数c=4*quadl(f,0,pi/2); %用数值积分求周长c1=[c1,c]; %连接周长end %结束循环syms x k %定义符号变量s=sqrt(1-k*sin(x)^2); %被积符号函数i=4*int(s,0,pi/2); %符号积分c2=subs(i,k,e.^2); %周长phi=linspace(0,pi/2); %角度向量[E,PHI]=meshgrid(e,phi); %矩阵F=4*sqrt(1-E.^2.*cos(PHI).^2); %被积函数c3=trapz(F)*phi(2); %用梯形法求周长[K,E]=ellipke(e.^2); %两类完全椭圆积分c4=4*E; %周长fs=16; %字体大小figure %创建图形窗口plot(e,c1,e,c3,'.',e,c3,'o',e,c4,'h','MarkerSize',16)%画周长曲线grid on%加网格title('椭圆的周长','fontsize',fs) %显示标题xlabel('\ite','fontsize',fs) %显示横坐标ylabel('\itC/a','fontsize',fs) %显示纵坐标legend('数值积分','符号积分','梯形求和','完全椭圆积分')%图例[图示]如zqy4.3.1图所示，椭圆的周长随偏心率的增加而单调减小。

当e= 0时，椭圆演化成没有偏心率的圆，周长为C = 2πa。

当e = 1时，椭圆退化成没有面积的线段，周长为C = 4a。

zqy4.3.1图zqy4.3.2图将周长的椭圆积分公式按二项式展开得π/221(21)!!4[1(sin)]d2!(21)nnnnC a en nϕϕ∞=-=--∑⎰其中(2n– 1)!! = 1·3…·(2n– 1)。

利用积分π/22(21)!!πsin d2!2nnnx xn-=⎰可得用无穷级数表示的椭圆周长211(21)!!2π{1[]}212!nnnnC a en n∞=-=--∑[程序]zqy4_3_2.m如下。

eps=1e-4; %容差nm=0; %圆的周长的项数为零c5=2*pi; %圆的周长for i=2:length(e) %按偏心率循环s=2*pi; %常数项之值for n=1:150 %按项数循环sn=-2*pi*(prod(1:2:2*n-1)*e(i)^n/2^n/factorial(n))^2/(2*n-1);%求各项之值 s=s+sn; %累加各项之值if abs(s-c4(i))<eps,break,end%与精确值相差很小时退出循环end%结束循环c5=[c5,s]; %连接周长nm=[nm,n]; %连接项数end%结束循环figure %创建图形窗口plot(e,c4,e,c5,'r.') %画周长曲线grid on%加网格title('椭圆的周长','fontsize',fs) %显示标题xlabel('\ite','fontsize',fs) %显示横坐标ylabel('\itC/a','fontsize',fs) %显示纵坐标legend('完全椭圆积分','级数求和') %图例text(e,c5,num2str(nm'),'fontsize',fs) %显示项数text(0,5,['\it\epsilon\rm=',num2str(eps)],'fontsize',fs)%显示精度作业：求余弦函数的积分。