第四章 微积分问题的计算机求解
? 微积分问题的解析解
? 函数的级数展开与级数求和问题求解 *
? 数值微分
? 数值积分问题
? 曲线积分与曲面积分的计算 *
4.1 微积分问题的解析解
4.1.1 极限问题的解析解
? 单变量函数的极限
– 格式 1,L= limit( fun,x,x0)
– 格式 2,L= limit( fun,x,x0,‘left’ 或 ‘ right’)
? 例, 试求解极限问题
>> syms x a b;
>> f=x*(1+a/x)^x*sin(b/x);
>> L=limit(f,x,inf)
L =
exp(a)*b
? 例,求解单边极限问题
>> syms x;
>> limit((exp(x^3)-1)/(1-cos(sqrt(x-sin(x)))),x,0,'right')
ans =
12
? 在 (-0.1,0.1)区间绘制出函数曲线,
>> x=-0.1:0.001:0.1;
>> y=(exp(x.^3)-1)./(1-cos(sqrt(x-sin(x))));
Warning,Divide by zero,
(Type "warning off
MATLAB,
divideByZero" to
suppress this warning.)
>> plot(x,y,'-',[0],
[12],'o')
? 多变量函数的极限,
– 格式,L1=limit(limit(f,x,x0),y,y0)
或 L1=limit(limit(f,y,y0),x,x0)
如果 x0 或 y0不是确定的值，而是另一个
变量的函数，如 x->g(y)，则上述的极限求
取顺序不能交换。
? 例：求出二元函数极限值
>> syms x y a;
>> f=exp(-
1/(y^2+x^2))*sin(x)^2/x^2*(1+1/y^2)^(x+a^2*y^2);
>> L=limit(limit(f,x,1/sqrt(y)),y,inf)
L =
exp(a^2)
4.1.2 函数导数的解析解
? 函数的导数和高阶导数
– 格式,y=diff(fun,x) %求导数
y= diff(fun,x,n) %求 n阶导数
? 例,
一阶导数,
>> syms x; f=sin(x)/(x^2+4*x+3);
>> f1=diff(f); pretty(f1)
cos(x) sin(x) (2 x + 4)
--------------- - -------------------
2 2 2
x + 4 x + 3 (x + 4 x + 3)
原函数及一阶导数图,
>> x1=0:.01:5;
>> y=subs(f,x,x1);
>> y1=subs(f1,x,x1);
>> plot(x1,y,x1,y1,‘:’)
更高阶导数,
>> tic,diff(f,x,100); toc
elapsed_time =
4.6860
? 原函数 4阶导数
>> f4=diff(f,x,4); pretty(f4)
2 sin(x) cos(x) (2 x + 4) sin(x) (2 x + 4)
------------ + 4 ------------------- - 12 ----------------- 2 2 2 2 3
x + 4 x + 3 (x + 4 x + 3) (x + 4 x + 3)
3 sin(x) cos(x) (2 x + 4) cos(x) (2 x + 4)
+ 12 --------------- - 24 ----------------- + 48 ---------------- 2 2 2 4 2 3
(x + 4 x + 3) (x + 4 x + 3) (x + 4 x + 3)
4 2 sin(x) (2 x + 4) sin(x) (2 x + 4) sin(x)
+ 24 ----------------- - 72 ----------------- + 24 --------------- 2 5 2 4 2 3
(x + 4 x + 3) (x + 4 x + 3) (x + 4 x + 3)
? 多元函数的偏导,
– 格式,f=diff(diff(f,x,m),y,n)
或 f=diff(diff(f,y,n),x,m)
? 例,求其偏导数并
用图表示。
>> syms x y z=(x^2-2*x)*exp(-x^2-y^2-x*y);
>> zx=simple(diff(z,x))
zx =
-exp(-x^2-y^2-x*y)*(-2*x+2+2*x^3+x^2*y-4*x^2-2*x*y)
>> zy=diff(z,y)
zy =
(x^2-2*x)*(-2*y-x)*exp(-x^2-y^2-x*y)
? 直接绘制三维曲面
>> [x,y]=meshgrid(-3:.2:3,-2:.2:2);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,z),axis([-3 3 -2 2 -0.7 1.5])
>> contour(x,y,z,30),hold on % 绘制等值线
>> zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-
4*x.^2-2*x.*y);
>> zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y); % 偏导
的数值解
>> quiver(x,y,zx,zy) % 绘制引力线
? 例
>> syms x y z; f=sin(x^2*y)*exp(-x^2*y-z^2);
>> df=diff(diff(diff(f,x,2),y),z); df=simple(df);
>> pretty(df)
2 2 2 2 2
-4 z exp(-x y - z ) (cos(x y) - 10 cos(x y) y x + 4
2 4 2 2 4 2 2
sin(x y) x y+ 4 cos(x y) x y - sin(x y))
? 多元函数的 Jacobi矩阵,
– 格式,J=jacobian(Y,X)
其中,X是自变量构成的向量,Y是由各个函数构成的
向量。
? 例,
试推导其 Jacobi 矩阵
>> syms r theta phi;
>> x=r*sin(theta)*cos(phi);
>> y=r*sin(theta)*sin(phi);
>> z=r*cos(theta);
>> J=jacobian([x; y; z],[r theta phi])
J =
[ sin(theta)*cos(phi),r*cos(theta)*cos(phi),-r*sin(theta)*sin(phi)]
[ sin(theta)*sin(phi),r*cos(theta)*sin(phi),r*sin(theta)*cos(phi)]
[ cos(theta),-r*sin(theta),0]
? 隐函数的偏导数,
– 格式,F=-diff(f,xj)/diff(f,xi)
? 例,
>> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y);
>> pretty(-simple(diff(f,x)/diff(f,y)))
3 2 2
-2 x + 2 + 2 x + x y - 4 x - 2 x y
- -----------------------------------------
x (x - 2) (2 y + x)
? 参数方程的导数
– 已知参数方程,求
– 格式,diff(f,t,k)/diff(g,t,k)
? 例,
>> syms t; y=sin(t)/(t+1)^3; x=cos(t)/(t+1)^3;
>> pretty(diff(y,t,4)/diff(x,t,4))
4.1.3 积分问题的解析解
? 不定积分的推导,
– 格式,F=int(fun,x)
? 例,
用 diff() 函数求其一阶导数，再积分,检验是否可以
得出一致的结果。
>> syms x; y=sin(x)/(x^2+4*x+3); y1=diff(y);
>> y0=int(y1); pretty(y0) % 对导数积分
sin(x) sin(x)
- 1/2 ------ + 1/2 ------
x + 3 x + 1
? 对原函数求 4 阶导数，再对结果进行 4次积分
>> y4=diff(y,4);
>> y0=int(int(int(int(y4))));
>> pretty(simple(y0))
sin(x)
------------
2
x + 4 x + 3
? 例,证明
>> syms a x; f=simple(int(x^3*cos(a*x)^2,x))
f = 1/16*(4*a^3*x^3*sin(2*a*x)+2*a^4
*x^4+6*a^2*x^2*cos(2*a*x)-6*a*x*sin(2*a*x)-
3*cos(2*a*x)-3)/a^4
>> f1=x^4/8+(x^3/(4*a)-3*x/(8*a^3))*sin(2*a*x)+..,
(3*x^2/(8*a^2)-3/(16*a^4))*cos(2*a*x);
>> simple(f-f1) % 求两个结果的差
ans =
-3/16/a^4
? 定积分与无穷积分计算,
– 格式,I=int(f,x,a,b)
– 格式,I=int(f,x,a,inf)
? 多重积分问题的 MATLAB求解
? 例,
>> syms x y z; f0=-4*z*exp(-x^2*y-z^2)*(cos(x^2*y)-
10*cos(x^2*y)*y*x^2+..,
4*sin(x^2*y)*x^4*y^2+4*cos(x^2*y)*x^4*y^2-sin(x^2*y));
>> f1=int(f0,z);f1=int(f1,y);f1=int(f1,x);
>> f1=simple(int(f1,x))
f1 =
exp(-x^2*y-z^2)*sin(x^2*y)
>> f2=int(f0,z); f2=int(f2,x); f2=int(f2,x);
>> f2=simple(int(f2,y))
f2 =
2*exp(-x^2*y-z^2)*tan(1/2*x^2*y)/(1+tan(1/2*x^2*y)^2)
>> simple(f1-f2)
ans =
0
顺序的改变使化简结果不同于原函数，但
其误差为 0，表明二者实际完全一致。这是由
于积分顺序不同，得不出实际的最简形式。
? 例,
>> syms x y z
>> int(int(int(4*x*z*exp(-x^2*y-z^2),x,0,1),y,0,pi),z,0,pi)
ans =
(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2))*pi^2*hypergeo
m([1],[2],-pi^2)
Ei(n,z)为指数积分，无解析解，但可求其数值解,
>> vpa(ans,60)
ans =
3.10807940208541272283461464767138521019142306317
021863483588
4.2 数值微分
4.2.1 数值微分算法
两种中心差分,
4.2.2 中心差分方法及其 MATLAB 实现
function [dy,dx]=diff_ctr(y,Dt,n)
yx1=[y 0 0 0 0 0]; yx2=[0 y 0 0 0 0]; yx3=[0 0 y 0 0 0];
yx4=[0 0 0 y 0 0]; yx5=[0 0 0 0 y 0]; yx6=[0 0 0 0 0 y];
switch n
case 1
dy = (-diff(yx1)+7*diff(yx2)+7*diff(yx3)- …
diff(yx4))/(12*Dt); L0=3;
case 2
dy=(-diff(yx1)+15*diff(yx2)- 15*diff(yx3)…
+diff(yx4))/(12*Dt^2);L0=3;
case 3
dy=(-diff(yx1)+7*diff(yx2)-6*diff(yx3)-6*diff(yx4)+..,
7*diff(yx5)-diff(yx6))/(8*Dt^3); L0=5;
case 4
dy = (-diff(yx1)+11*diff(yx2)-28*diff(yx3)+28*…
diff(yx4)-11*diff(yx5)+diff(yx6))/(6*Dt^4); L0=5;
end
dy=dy(L0+1:end-L0); dx=([1:length(dy)]+L0-2-(n>2))*Dt;
调用格式,
y为 等距实测数据,dy为得出的导数向量,dx为
相应的自变量向量 。
[,] _ (,,)yxd d d i f f c t r y t n??
? 例,
求导数的解析解,再 用数值微分求取原函数的 1～
4 阶导数，并和解析解比较精度。
>> h=0.05; x=0:h:pi;
>> syms x1; y=sin(x1)/(x1^2+4*x1+3);
% 求各阶导数的解析解与对照数据
>> yy1=diff(y); f1=subs(yy1,x1,x);
>> yy2=diff(yy1); f2=subs(yy2,x1,x);
>> yy3=diff(yy2); f3=subs(yy3,x1,x);
>> yy4=diff(yy3); f4=subs(yy4,x1,x);
>> y=sin(x)./(x.^2+4*x+3); % 生成已知数据点
>> [y1,dx1]=diff_ctr(y,h,1); subplot(221),plot(x,f1,dx1,y1,':');
>> [y2,dx2]=diff_ctr(y,h,2); subplot(222),plot(x,f2,dx2,y2,':')
>> [y3,dx3]=diff_ctr(y,h,3); subplot(223),plot(x,f3,dx3,y3,':');
>> [y4,dx4]=diff_ctr(y,h,4); subplot(224),plot(x,f4,dx4,y4,':')
求最大相对误差,
>> norm((y4-…
f4(4:60))./f4(4:60))
ans =
3.5025e-004
4.2.3 插值多项式的导数
? 基本思想：当已知函数在一些离散点上的函
数值时，该函数可用插值多项式来近似，然
后对多项式进行微分求得导数。
? 选取 x=0附近的少量点
? 进行插值多项式拟合
? g(x)在 x=0处的 k阶导数为
(,),1,2,,1iix y i n??L
1
1 2 1()
nn
nng x c x c x c x c
?
?? ? ? ? ?L
() 1( 0 ) ! 0,1,2,,k nkg c k k n???? L
? 通过坐标变换用上述方法计算任意 x点处的导
数值
? 令
? 将 g(x)写成 z的表达式
? 导数为
? 可直接用 拟合节点 得到系数
d=polyfit(x-a,y,length(xd)-1)
z x a??
1
1 2 1( ) ( )
nn
nng x g z d z d z d z d
?
?? ? ? ? ? ?L
( ) ( )
1( ) ( 0 ) ! 0,1,,
kk
nkg a g d k k n??? ? ? L
()gz (,)
iix a y? i
d
? 例：数据集合如下,
xd,0 0.2000 0.4000 0.6000 0.8000 1.000
yd,0.3927 0.5672 0.6982 0.7941 0.8614 0.9053
计算 x=a=0.3处的各阶导数。
>> xd=[ 0 0.2000 0.4000 0.6000 0.8000 1.000];
>> yd=[0.3927 0.5672 0.6982 0.7941 0.8614 0.9053];
>> a=0.3;L=length(xd);
>> d=polyfit(xd-a,yd,L-1);fact=[1];
>> for k=1:L-1;fact=[factorial(k),fact];end
>> deriv=d.*fact
deriv =
1.8750 -1.3750 1.0406 -0.9710 0.6533 0.6376
? 建立计算插值多项式各阶导数的 poly_drv.m
function der=poly_drv(xd,yd,a)
m=length(xd)-1;
d=polyfit(xd-a,yd,m);
c=d(m:-1:1);
fact(1)=1;for i=2:m; fact(i)=i*fact(i-1);end
der=c.*fact;
? 例,
>> xd=[ 0 0.2000 0.4000 0.6000 0.8000 1.000];
>> yd=[0.3927 0.5672 0.6982 0.7941 0.8614 0.9053];
>> a=0.3; der=poly_drv(xd,yd,a)
der =
0.6533 -0.9710 1.0406 -1.3750 1.8750
4.2.4 二元函数的梯度计算
? 格式,
? 若 z矩阵是建立在等间距的形式生成的网格基
础上，则实际梯度为
[,] ( )xyf f g r a d i e n t z?
/,/x x y yf f x f f y? ? ? ?
? 例,
计算梯度，绘制引力线图,
>> [x,y]=meshgrid(-3:.2:3,-2:.2:2); z=(x.^2-2*x).*exp(-
x.^2-y.^2-x.*y);
>> [fx,fy]=gradient(z);
>> fx=fx/0.2; fy=fy/0.2;
>> contour(x,y,z,30);
>> hold on;
>> quiver(x,y,fx,fy)
%绘制等高线与
引力线图
? 绘制误差曲面,
>> zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-
4*x.^2-2*x.*y);
>> zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,abs(fx-zx)); axis([-3 3 -2 2 0,0.08])
>> figure; surf(x,y,abs(fy-zy)); axis([-3 3 -2 2 0,0.11])
? 为减少误差，对网格加密一倍,
>> [x,y]=meshgrid(-3:.1:3,-2:.1:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> [fx,fy]=gradient(z); fx=fx/0.1; fy=fy/0.1;
>> zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-4*x.^2-2*x.*y);
>> zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,abs(fx-zx)); axis([-3 3 -2 2 0,0.02])
>> figure; surf(x,y,abs(fy-zy)); axis([-3 3 -2 2 0,0.06])
4.3 数值积分问题
4.3.1 由给定数据进行梯形求积
Sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2
? 格式,S=trapz(x,y)
? 例,
>> x1=[0:pi/30:pi]'; y=[sin(x1) cos(x1) sin(x1/2)];
>> x=[x1 x1 x1]; S=sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2
S =
1.9982 0.0000 1.9995
>> S1=trapz(x1,y) % 得出和上述完全一致的结果
S1 =
1.9982 0.0000 1.9995
? 例,
画图
>> x=[0:0.01:3*pi/2,3*pi/2]; % 这样赋值能确保 3*pi/2点
被包含在内
>> y=cos(15*x); plot(x,y)
% 求取理论值
>> syms x,A=int(cos…
(15*x),0,3*pi/2)
A =
1/15
随着步距 h的减小，计算精度逐渐增加,
>> h0=[0.1,0.01,0.001,0.0001,0.00001,0.000001]; v=[];
>> for h=h0,
x=[0:h:3*pi/2,3*pi/2]; y=cos(15*x); I=trapz(x,y);
v=[v; h,I,1/15-I ];
end
>> v
v =
0.1000 0.0539 0.0128
0.0100 0.0665 0.0001
0.0010 0.0667 0.0000
0.0001 0.0667 0.0000
0.0000 0.0667 0.0000
0.0000 0.0667 0.0000
4.3.2 单变量数值积分问题求解
? 格式,
y=quad(Fun,a,b) y=quadl(Fun,a,b) % 求定积分
y=quad(Fun,a,b,) y=quadl(Fun,a,b,) %限定精
度的定积分求解，默认精度为 10－ 6。
? ?
? 例,
第三种：匿名函数 (MATLAB 7.0)
第二种,inline 函数
第一种，一般函数方法
函数定义被积函数,
>> y=quad('c3ffun',0,1.5)
y =
0.9661
? 用 inline 函数定义被积函数,
>> f=inline('2/sqrt(pi)*exp(-x.^2)','x');
>> y=quad(f,0,1.5)
y =
0.9661
? 运用符号工具箱,
>> syms x,y0=vpa(int(2/sqrt(pi)*exp(-x^2),0,1.5),60)
y0 =
,96610514647531071393693372994990579499622494325746147328
5749
>> y=quad(f,0,1.5,1e-20) % 设置高精度，但该方法失效
? 例,
提高求解精度,
>> y=quadl(f,0,1.5,1e-20)
y =
0.9661
>> abs(y-y0)
ans =
,6402522848913892e-16
>> format long ％ 16位精度
>> y=quadl(f,0,1.5,1e-20)
y =
0.96610514647531
? 例,求解
绘制函数,
>> x=[0:0.01:2,2+eps:0.01:4,4];
>> y=exp(x.^2).*(x<=2)+80./(4-sin(16*pi*x)).*(x>2);
>> y(end)=0;
>> x=[eps,x];
>> y=[0,y];
>> fill(x,y,'g')
? 调用 quad( ),
>> f=inline('exp(x.^2).*(x<=2)+80*(x>2)./(4-sin(16*pi*x))','x');
>> I1=quad(f,0,4)
I1 =
57.76435412500863
? 调用 quadl( ),
>> I2=quadl(f,0,4)
I2 =
57.76445016946768
>> syms x; I=vpa(int(exp(x^2),0,2)+int(80/(4-sin(16*pi*x)),2,4))
I =
57.764450125053010333315235385182
4.3.3 Gauss求积公式
? 为使求积公式得到较高的代数精度
? 对求积区间 [a,b]，通过变换
? 有
1
1 0
( ) ( )
n
kk
k
f x d x A f x
? ?
? ??
22
b a b axt????
1
1 0
( ) ( ) ( )2 2 2 2 2 2
nb
ka
k
b a b a a b b a b a a bf x d x f t d t A f t
? ?
? ? ? ? ? ?? ? ? ????
1,0, 5 7 7 3 5 0 3,1, 0 0 0 0 0 0 0 0 ;
2,0, 7 7 4 5 9 6 7,0, 5 5 5 5 5 5 5 6
0, 0 0 0 0 0 0 0 0,0, 8 8 8 8 8 8 8 9 ;
kk
kk
kk
n x A
n x A
xA
? ? ? ?
? ? ? ?
??
? 以 n=2的高斯公式为例,
function g=gauss2(fun,a,b)
h=(b-a)/2;
c=(a+b)/2;
x=[h*(-0.7745967)+c,c,h*0.7745967+c];
g=h*(0.55555556*(gaussf(x(1))+gaussf(x(3)))+0.8888888
9*gaussf(x(2)));
function y=gaussf(x)
y=cos(x);
>> gauss2('gaussf',0,1)
ans =
0.8415
4.3.4 基于样条插值的数值微积分运算
? 基于样条插值的数值微分运算
– 格式,
Sd=fnder(S,k)
该函数可以求取 S的 k阶导数。
– 格式,
Sd=fnder(S,[k1,…,kn])
可以求取多变量函数的偏导数
? 例,
>> syms x; f=(x^2-3*x+5)*exp(-5*x)*sin(x);
>> ezplot(diff(f),[0,1]),hold on
>> x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> sp1=csapi(x,y);
>> dsp1=fnder(sp1,1);
>> fnplt(dsp1,'--')
>> sp2=spapi(5,x,y);
>> dsp2=fnder(sp2,1);
>> fnplt(dsp2,':');
>> axis([0,1,-0.8,5])
? 例,
拟合曲面
>> x0=-3:.3:3; y0=-2:.2:2; [x,y]=ndgrid(x0,y0);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> sp=spapi({5,5},…
{x0,y0},z);
>>dspxy=fnder(sp,[1,1]);
>> fnplt(dspxy)
理论方法,
>> syms x y; z=(x^2-2*x)*exp(-x^2-y^2-x*y);
>> ezsurf(diff(diff(z,x),y),[-3 3],[-2 2])
? 基于样条插值的数值积分运算
– 格式,
f=fnint(S)
其中 S为样条函数。
例,考虑 中较稀疏的样本点,用样
条积分的方式求出定积分及积分函数。
>> x=[0,0.4,1 2,pi]; y=sin(x);
>> sp1=csapi(x,y); a=fnint(sp1,1);
>> xx=fnval(a,[0,pi]); xx(2)-xx(1)
ans =
2.0191
0 s in xd x
??
>> sp2=spapi(5,x,y); b=fnint(sp2,1);
>> xx=fnval(b,[0,pi]); xx(2)-xx(1)
ans =
1.9999
绘制曲线
>> ezplot('-cos(t)+2',[0,pi]); hold on
>> fnplt(a,'--');
>> fnplt(b,':')
4.3.5 双重积分问题的数值解
? 矩形区域上的二重积分的数值计算
? 格式,
矩形区域的双重积分,
y=dblquad(Fun,xm,xM,ym,yM)
限定精度的双重积分,
y=dblquad(Fun,xm,xM,ym,yM,) ?
? 例,求解
>> f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');
>> y=dblquad(f,-2,2,-1,1)
y =
1.57449318974494
? 任意区域上二元函数的数值积分 (调用工具箱）
? 格式,
一般双重积分
J= quad2dggen(fun,ymin,ymax,xlower,xupper)
限定精度的双重积分
J= quad2dggen(fun,ymin,ymax,xlower,xupper,)
?
例
>> fh=inline('sqrt(1-x.^2/2)','x'); % 内积分上限
>> fl=inline('-sqrt(1-x.^2/2)','x'); % 内积分下限
>> f=inline('exp(-x.^2/2).*sin(x.^2+y)','y','x'); %
交换顺序的被积函数
>> y=quad2dggen(f,fl,fh,-1/2,1,eps)
y =
0.41192954617630
? 解析解方法,
>> syms x y
>> i1=int(exp(-x^2/2)*sin(x^2+y),y,-sqrt(1-x^2/2),sqrt(1-
x^2/2));
>> int(i1,x,-1/2,1)
Warning,Explicit integral could not be found,
> In D:\MATLAB6p5\toolbox\symbolic\@sym\int.m at line
58
ans =
int(2*exp(-1/2*x^2)*sin(x^2)*sin(1/2*(4-2*x^2)^(1/2)),x =
-1/2,,1)
>> vpa(ans)
ans =
.41192954617629511965175994017601
2 2
2
11
22
1 1
s i n ( )
y x
y
I e x y d x d y
?
?
? ??
????
2
22
22
1
sin ( )
x
xy
I e x y d x d y
?
??
????
例,计算单位圆域上的积分,
先把二重积分转化为二次积分的形式,
>> syms x y
i1=int(exp(-x^2/2)*sin(x^2+y),x,-sqrt(1-y.^2),sqrt(1-y.^2));
Warning,Explicit integral could not be found,
> In D:\MATLAB6p5\toolbox\symbolic\@sym\int.m at line 58
对 x是不可积的，故调用解析解方法不会得出结
果，而数值解求解不受此影响。
>> fh=inline('sqrt(1-y.^2)','y'); % 内积分上限
>> fl=inline('-sqrt(1-y.^2)','y'); % 内积分下限
>> f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y'); %交换顺序
的被积函数
>> I=quad2dggen(f,fl,fh,-1,1,eps)
Integral did not converge--singularity likely
I =
0.53686038269795
4.3.6 三重定积分的数值求解
? 格式,
I=triplequad(Fun,xm,xM,ym,yM,zm,zM，,@quadl)
其中 @quadl为具体求解一元积分的数值函数,
也可选用 @quad或自编积分函数,但格式要与
quadl一致。
?
? 例,
>> triplequad(inline('4*x.*z.*exp(-x.*x.*y-z.*z)',…
'x','y','z'),0,2,0,pi,0,pi,1e-7,@quadl)
ans =
3.10807945143834
? 微积分问题的解析解
? 函数的级数展开与级数求和问题求解 *
? 数值微分
? 数值积分问题
? 曲线积分与曲面积分的计算 *
4.1 微积分问题的解析解
4.1.1 极限问题的解析解
? 单变量函数的极限
– 格式 1,L= limit( fun,x,x0)
– 格式 2,L= limit( fun,x,x0,‘left’ 或 ‘ right’)
? 例, 试求解极限问题
>> syms x a b;
>> f=x*(1+a/x)^x*sin(b/x);
>> L=limit(f,x,inf)
L =
exp(a)*b
? 例,求解单边极限问题
>> syms x;
>> limit((exp(x^3)-1)/(1-cos(sqrt(x-sin(x)))),x,0,'right')
ans =
12
? 在 (-0.1,0.1)区间绘制出函数曲线,
>> x=-0.1:0.001:0.1;
>> y=(exp(x.^3)-1)./(1-cos(sqrt(x-sin(x))));
Warning,Divide by zero,
(Type "warning off
MATLAB,
divideByZero" to
suppress this warning.)
>> plot(x,y,'-',[0],
[12],'o')
? 多变量函数的极限,
– 格式,L1=limit(limit(f,x,x0),y,y0)
或 L1=limit(limit(f,y,y0),x,x0)
如果 x0 或 y0不是确定的值，而是另一个
变量的函数，如 x->g(y)，则上述的极限求
取顺序不能交换。
? 例：求出二元函数极限值
>> syms x y a;
>> f=exp(-
1/(y^2+x^2))*sin(x)^2/x^2*(1+1/y^2)^(x+a^2*y^2);
>> L=limit(limit(f,x,1/sqrt(y)),y,inf)
L =
exp(a^2)
4.1.2 函数导数的解析解
? 函数的导数和高阶导数
– 格式,y=diff(fun,x) %求导数
y= diff(fun,x,n) %求 n阶导数
? 例,
一阶导数,
>> syms x; f=sin(x)/(x^2+4*x+3);
>> f1=diff(f); pretty(f1)
cos(x) sin(x) (2 x + 4)
--------------- - -------------------
2 2 2
x + 4 x + 3 (x + 4 x + 3)
原函数及一阶导数图,
>> x1=0:.01:5;
>> y=subs(f,x,x1);
>> y1=subs(f1,x,x1);
>> plot(x1,y,x1,y1,‘:’)
更高阶导数,
>> tic,diff(f,x,100); toc
elapsed_time =
4.6860
? 原函数 4阶导数
>> f4=diff(f,x,4); pretty(f4)
2 sin(x) cos(x) (2 x + 4) sin(x) (2 x + 4)
------------ + 4 ------------------- - 12 ----------------- 2 2 2 2 3
x + 4 x + 3 (x + 4 x + 3) (x + 4 x + 3)
3 sin(x) cos(x) (2 x + 4) cos(x) (2 x + 4)
+ 12 --------------- - 24 ----------------- + 48 ---------------- 2 2 2 4 2 3
(x + 4 x + 3) (x + 4 x + 3) (x + 4 x + 3)
4 2 sin(x) (2 x + 4) sin(x) (2 x + 4) sin(x)
+ 24 ----------------- - 72 ----------------- + 24 --------------- 2 5 2 4 2 3
(x + 4 x + 3) (x + 4 x + 3) (x + 4 x + 3)
? 多元函数的偏导,
– 格式,f=diff(diff(f,x,m),y,n)
或 f=diff(diff(f,y,n),x,m)
? 例,求其偏导数并
用图表示。
>> syms x y z=(x^2-2*x)*exp(-x^2-y^2-x*y);
>> zx=simple(diff(z,x))
zx =
-exp(-x^2-y^2-x*y)*(-2*x+2+2*x^3+x^2*y-4*x^2-2*x*y)
>> zy=diff(z,y)
zy =
(x^2-2*x)*(-2*y-x)*exp(-x^2-y^2-x*y)
? 直接绘制三维曲面
>> [x,y]=meshgrid(-3:.2:3,-2:.2:2);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,z),axis([-3 3 -2 2 -0.7 1.5])
>> contour(x,y,z,30),hold on % 绘制等值线
>> zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-
4*x.^2-2*x.*y);
>> zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y); % 偏导
的数值解
>> quiver(x,y,zx,zy) % 绘制引力线
? 例
>> syms x y z; f=sin(x^2*y)*exp(-x^2*y-z^2);
>> df=diff(diff(diff(f,x,2),y),z); df=simple(df);
>> pretty(df)
2 2 2 2 2
-4 z exp(-x y - z ) (cos(x y) - 10 cos(x y) y x + 4
2 4 2 2 4 2 2
sin(x y) x y+ 4 cos(x y) x y - sin(x y))
? 多元函数的 Jacobi矩阵,
– 格式,J=jacobian(Y,X)
其中,X是自变量构成的向量,Y是由各个函数构成的
向量。
? 例,
试推导其 Jacobi 矩阵
>> syms r theta phi;
>> x=r*sin(theta)*cos(phi);
>> y=r*sin(theta)*sin(phi);
>> z=r*cos(theta);
>> J=jacobian([x; y; z],[r theta phi])
J =
[ sin(theta)*cos(phi),r*cos(theta)*cos(phi),-r*sin(theta)*sin(phi)]
[ sin(theta)*sin(phi),r*cos(theta)*sin(phi),r*sin(theta)*cos(phi)]
[ cos(theta),-r*sin(theta),0]
? 隐函数的偏导数,
– 格式,F=-diff(f,xj)/diff(f,xi)
? 例,
>> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y);
>> pretty(-simple(diff(f,x)/diff(f,y)))
3 2 2
-2 x + 2 + 2 x + x y - 4 x - 2 x y
- -----------------------------------------
x (x - 2) (2 y + x)
? 参数方程的导数
– 已知参数方程,求
– 格式,diff(f,t,k)/diff(g,t,k)
? 例,
>> syms t; y=sin(t)/(t+1)^3; x=cos(t)/(t+1)^3;
>> pretty(diff(y,t,4)/diff(x,t,4))
4.1.3 积分问题的解析解
? 不定积分的推导,
– 格式,F=int(fun,x)
? 例,
用 diff() 函数求其一阶导数，再积分,检验是否可以
得出一致的结果。
>> syms x; y=sin(x)/(x^2+4*x+3); y1=diff(y);
>> y0=int(y1); pretty(y0) % 对导数积分
sin(x) sin(x)
- 1/2 ------ + 1/2 ------
x + 3 x + 1
? 对原函数求 4 阶导数，再对结果进行 4次积分
>> y4=diff(y,4);
>> y0=int(int(int(int(y4))));
>> pretty(simple(y0))
sin(x)
------------
2
x + 4 x + 3
? 例,证明
>> syms a x; f=simple(int(x^3*cos(a*x)^2,x))
f = 1/16*(4*a^3*x^3*sin(2*a*x)+2*a^4
*x^4+6*a^2*x^2*cos(2*a*x)-6*a*x*sin(2*a*x)-
3*cos(2*a*x)-3)/a^4
>> f1=x^4/8+(x^3/(4*a)-3*x/(8*a^3))*sin(2*a*x)+..,
(3*x^2/(8*a^2)-3/(16*a^4))*cos(2*a*x);
>> simple(f-f1) % 求两个结果的差
ans =
-3/16/a^4
? 定积分与无穷积分计算,
– 格式,I=int(f,x,a,b)
– 格式,I=int(f,x,a,inf)
? 多重积分问题的 MATLAB求解
? 例,
>> syms x y z; f0=-4*z*exp(-x^2*y-z^2)*(cos(x^2*y)-
10*cos(x^2*y)*y*x^2+..,
4*sin(x^2*y)*x^4*y^2+4*cos(x^2*y)*x^4*y^2-sin(x^2*y));
>> f1=int(f0,z);f1=int(f1,y);f1=int(f1,x);
>> f1=simple(int(f1,x))
f1 =
exp(-x^2*y-z^2)*sin(x^2*y)
>> f2=int(f0,z); f2=int(f2,x); f2=int(f2,x);
>> f2=simple(int(f2,y))
f2 =
2*exp(-x^2*y-z^2)*tan(1/2*x^2*y)/(1+tan(1/2*x^2*y)^2)
>> simple(f1-f2)
ans =
0
顺序的改变使化简结果不同于原函数，但
其误差为 0，表明二者实际完全一致。这是由
于积分顺序不同，得不出实际的最简形式。
? 例,
>> syms x y z
>> int(int(int(4*x*z*exp(-x^2*y-z^2),x,0,1),y,0,pi),z,0,pi)
ans =
(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2))*pi^2*hypergeo
m([1],[2],-pi^2)
Ei(n,z)为指数积分，无解析解，但可求其数值解,
>> vpa(ans,60)
ans =
3.10807940208541272283461464767138521019142306317
021863483588
4.2 数值微分
4.2.1 数值微分算法
两种中心差分,
4.2.2 中心差分方法及其 MATLAB 实现
function [dy,dx]=diff_ctr(y,Dt,n)
yx1=[y 0 0 0 0 0]; yx2=[0 y 0 0 0 0]; yx3=[0 0 y 0 0 0];
yx4=[0 0 0 y 0 0]; yx5=[0 0 0 0 y 0]; yx6=[0 0 0 0 0 y];
switch n
case 1
dy = (-diff(yx1)+7*diff(yx2)+7*diff(yx3)- …
diff(yx4))/(12*Dt); L0=3;
case 2
dy=(-diff(yx1)+15*diff(yx2)- 15*diff(yx3)…
+diff(yx4))/(12*Dt^2);L0=3;
case 3
dy=(-diff(yx1)+7*diff(yx2)-6*diff(yx3)-6*diff(yx4)+..,
7*diff(yx5)-diff(yx6))/(8*Dt^3); L0=5;
case 4
dy = (-diff(yx1)+11*diff(yx2)-28*diff(yx3)+28*…
diff(yx4)-11*diff(yx5)+diff(yx6))/(6*Dt^4); L0=5;
end
dy=dy(L0+1:end-L0); dx=([1:length(dy)]+L0-2-(n>2))*Dt;
调用格式,
y为 等距实测数据,dy为得出的导数向量,dx为
相应的自变量向量 。
[,] _ (,,)yxd d d i f f c t r y t n??
? 例,
求导数的解析解,再 用数值微分求取原函数的 1～
4 阶导数，并和解析解比较精度。
>> h=0.05; x=0:h:pi;
>> syms x1; y=sin(x1)/(x1^2+4*x1+3);
% 求各阶导数的解析解与对照数据
>> yy1=diff(y); f1=subs(yy1,x1,x);
>> yy2=diff(yy1); f2=subs(yy2,x1,x);
>> yy3=diff(yy2); f3=subs(yy3,x1,x);
>> yy4=diff(yy3); f4=subs(yy4,x1,x);
>> y=sin(x)./(x.^2+4*x+3); % 生成已知数据点
>> [y1,dx1]=diff_ctr(y,h,1); subplot(221),plot(x,f1,dx1,y1,':');
>> [y2,dx2]=diff_ctr(y,h,2); subplot(222),plot(x,f2,dx2,y2,':')
>> [y3,dx3]=diff_ctr(y,h,3); subplot(223),plot(x,f3,dx3,y3,':');
>> [y4,dx4]=diff_ctr(y,h,4); subplot(224),plot(x,f4,dx4,y4,':')
求最大相对误差,
>> norm((y4-…
f4(4:60))./f4(4:60))
ans =
3.5025e-004
4.2.3 插值多项式的导数
? 基本思想：当已知函数在一些离散点上的函
数值时，该函数可用插值多项式来近似，然
后对多项式进行微分求得导数。
? 选取 x=0附近的少量点
? 进行插值多项式拟合
? g(x)在 x=0处的 k阶导数为
(,),1,2,,1iix y i n??L
1
1 2 1()
nn
nng x c x c x c x c
?
?? ? ? ? ?L
() 1( 0 ) ! 0,1,2,,k nkg c k k n???? L
? 通过坐标变换用上述方法计算任意 x点处的导
数值
? 令
? 将 g(x)写成 z的表达式
? 导数为
? 可直接用 拟合节点 得到系数
d=polyfit(x-a,y,length(xd)-1)
z x a??
1
1 2 1( ) ( )
nn
nng x g z d z d z d z d
?
?? ? ? ? ? ?L
( ) ( )
1( ) ( 0 ) ! 0,1,,
kk
nkg a g d k k n??? ? ? L
()gz (,)
iix a y? i
d
? 例：数据集合如下,
xd,0 0.2000 0.4000 0.6000 0.8000 1.000
yd,0.3927 0.5672 0.6982 0.7941 0.8614 0.9053
计算 x=a=0.3处的各阶导数。
>> xd=[ 0 0.2000 0.4000 0.6000 0.8000 1.000];
>> yd=[0.3927 0.5672 0.6982 0.7941 0.8614 0.9053];
>> a=0.3;L=length(xd);
>> d=polyfit(xd-a,yd,L-1);fact=[1];
>> for k=1:L-1;fact=[factorial(k),fact];end
>> deriv=d.*fact
deriv =
1.8750 -1.3750 1.0406 -0.9710 0.6533 0.6376
? 建立计算插值多项式各阶导数的 poly_drv.m
function der=poly_drv(xd,yd,a)
m=length(xd)-1;
d=polyfit(xd-a,yd,m);
c=d(m:-1:1);
fact(1)=1;for i=2:m; fact(i)=i*fact(i-1);end
der=c.*fact;
? 例,
>> xd=[ 0 0.2000 0.4000 0.6000 0.8000 1.000];
>> yd=[0.3927 0.5672 0.6982 0.7941 0.8614 0.9053];
>> a=0.3; der=poly_drv(xd,yd,a)
der =
0.6533 -0.9710 1.0406 -1.3750 1.8750
4.2.4 二元函数的梯度计算
? 格式,
? 若 z矩阵是建立在等间距的形式生成的网格基
础上，则实际梯度为
[,] ( )xyf f g r a d i e n t z?
/,/x x y yf f x f f y? ? ? ?
? 例,
计算梯度，绘制引力线图,
>> [x,y]=meshgrid(-3:.2:3,-2:.2:2); z=(x.^2-2*x).*exp(-
x.^2-y.^2-x.*y);
>> [fx,fy]=gradient(z);
>> fx=fx/0.2; fy=fy/0.2;
>> contour(x,y,z,30);
>> hold on;
>> quiver(x,y,fx,fy)
%绘制等高线与
引力线图
? 绘制误差曲面,
>> zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-
4*x.^2-2*x.*y);
>> zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,abs(fx-zx)); axis([-3 3 -2 2 0,0.08])
>> figure; surf(x,y,abs(fy-zy)); axis([-3 3 -2 2 0,0.11])
? 为减少误差，对网格加密一倍,
>> [x,y]=meshgrid(-3:.1:3,-2:.1:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> [fx,fy]=gradient(z); fx=fx/0.1; fy=fy/0.1;
>> zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-4*x.^2-2*x.*y);
>> zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y);
>> surf(x,y,abs(fx-zx)); axis([-3 3 -2 2 0,0.02])
>> figure; surf(x,y,abs(fy-zy)); axis([-3 3 -2 2 0,0.06])
4.3 数值积分问题
4.3.1 由给定数据进行梯形求积
Sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2
? 格式,S=trapz(x,y)
? 例,
>> x1=[0:pi/30:pi]'; y=[sin(x1) cos(x1) sin(x1/2)];
>> x=[x1 x1 x1]; S=sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2
S =
1.9982 0.0000 1.9995
>> S1=trapz(x1,y) % 得出和上述完全一致的结果
S1 =
1.9982 0.0000 1.9995
? 例,
画图
>> x=[0:0.01:3*pi/2,3*pi/2]; % 这样赋值能确保 3*pi/2点
被包含在内
>> y=cos(15*x); plot(x,y)
% 求取理论值
>> syms x,A=int(cos…
(15*x),0,3*pi/2)
A =
1/15
随着步距 h的减小，计算精度逐渐增加,
>> h0=[0.1,0.01,0.001,0.0001,0.00001,0.000001]; v=[];
>> for h=h0,
x=[0:h:3*pi/2,3*pi/2]; y=cos(15*x); I=trapz(x,y);
v=[v; h,I,1/15-I ];
end
>> v
v =
0.1000 0.0539 0.0128
0.0100 0.0665 0.0001
0.0010 0.0667 0.0000
0.0001 0.0667 0.0000
0.0000 0.0667 0.0000
0.0000 0.0667 0.0000
4.3.2 单变量数值积分问题求解
? 格式,
y=quad(Fun,a,b) y=quadl(Fun,a,b) % 求定积分
y=quad(Fun,a,b,) y=quadl(Fun,a,b,) %限定精
度的定积分求解，默认精度为 10－ 6。
? ?
? 例,
第三种：匿名函数 (MATLAB 7.0)
第二种,inline 函数
第一种，一般函数方法
函数定义被积函数,
>> y=quad('c3ffun',0,1.5)
y =
0.9661
? 用 inline 函数定义被积函数,
>> f=inline('2/sqrt(pi)*exp(-x.^2)','x');
>> y=quad(f,0,1.5)
y =
0.9661
? 运用符号工具箱,
>> syms x,y0=vpa(int(2/sqrt(pi)*exp(-x^2),0,1.5),60)
y0 =
,96610514647531071393693372994990579499622494325746147328
5749
>> y=quad(f,0,1.5,1e-20) % 设置高精度，但该方法失效
? 例,
提高求解精度,
>> y=quadl(f,0,1.5,1e-20)
y =
0.9661
>> abs(y-y0)
ans =
,6402522848913892e-16
>> format long ％ 16位精度
>> y=quadl(f,0,1.5,1e-20)
y =
0.96610514647531
? 例,求解
绘制函数,
>> x=[0:0.01:2,2+eps:0.01:4,4];
>> y=exp(x.^2).*(x<=2)+80./(4-sin(16*pi*x)).*(x>2);
>> y(end)=0;
>> x=[eps,x];
>> y=[0,y];
>> fill(x,y,'g')
? 调用 quad( ),
>> f=inline('exp(x.^2).*(x<=2)+80*(x>2)./(4-sin(16*pi*x))','x');
>> I1=quad(f,0,4)
I1 =
57.76435412500863
? 调用 quadl( ),
>> I2=quadl(f,0,4)
I2 =
57.76445016946768
>> syms x; I=vpa(int(exp(x^2),0,2)+int(80/(4-sin(16*pi*x)),2,4))
I =
57.764450125053010333315235385182
4.3.3 Gauss求积公式
? 为使求积公式得到较高的代数精度
? 对求积区间 [a,b]，通过变换
? 有
1
1 0
( ) ( )
n
kk
k
f x d x A f x
? ?
? ??
22
b a b axt????
1
1 0
( ) ( ) ( )2 2 2 2 2 2
nb
ka
k
b a b a a b b a b a a bf x d x f t d t A f t
? ?
? ? ? ? ? ?? ? ? ????
1,0, 5 7 7 3 5 0 3,1, 0 0 0 0 0 0 0 0 ;
2,0, 7 7 4 5 9 6 7,0, 5 5 5 5 5 5 5 6
0, 0 0 0 0 0 0 0 0,0, 8 8 8 8 8 8 8 9 ;
kk
kk
kk
n x A
n x A
xA
? ? ? ?
? ? ? ?
??
? 以 n=2的高斯公式为例,
function g=gauss2(fun,a,b)
h=(b-a)/2;
c=(a+b)/2;
x=[h*(-0.7745967)+c,c,h*0.7745967+c];
g=h*(0.55555556*(gaussf(x(1))+gaussf(x(3)))+0.8888888
9*gaussf(x(2)));
function y=gaussf(x)
y=cos(x);
>> gauss2('gaussf',0,1)
ans =
0.8415
4.3.4 基于样条插值的数值微积分运算
? 基于样条插值的数值微分运算
– 格式,
Sd=fnder(S,k)
该函数可以求取 S的 k阶导数。
– 格式,
Sd=fnder(S,[k1,…,kn])
可以求取多变量函数的偏导数
? 例,
>> syms x; f=(x^2-3*x+5)*exp(-5*x)*sin(x);
>> ezplot(diff(f),[0,1]),hold on
>> x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> sp1=csapi(x,y);
>> dsp1=fnder(sp1,1);
>> fnplt(dsp1,'--')
>> sp2=spapi(5,x,y);
>> dsp2=fnder(sp2,1);
>> fnplt(dsp2,':');
>> axis([0,1,-0.8,5])
? 例,
拟合曲面
>> x0=-3:.3:3; y0=-2:.2:2; [x,y]=ndgrid(x0,y0);
>> z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
>> sp=spapi({5,5},…
{x0,y0},z);
>>dspxy=fnder(sp,[1,1]);
>> fnplt(dspxy)
理论方法,
>> syms x y; z=(x^2-2*x)*exp(-x^2-y^2-x*y);
>> ezsurf(diff(diff(z,x),y),[-3 3],[-2 2])
? 基于样条插值的数值积分运算
– 格式,
f=fnint(S)
其中 S为样条函数。
例,考虑 中较稀疏的样本点,用样
条积分的方式求出定积分及积分函数。
>> x=[0,0.4,1 2,pi]; y=sin(x);
>> sp1=csapi(x,y); a=fnint(sp1,1);
>> xx=fnval(a,[0,pi]); xx(2)-xx(1)
ans =
2.0191
0 s in xd x
??
>> sp2=spapi(5,x,y); b=fnint(sp2,1);
>> xx=fnval(b,[0,pi]); xx(2)-xx(1)
ans =
1.9999
绘制曲线
>> ezplot('-cos(t)+2',[0,pi]); hold on
>> fnplt(a,'--');
>> fnplt(b,':')
4.3.5 双重积分问题的数值解
? 矩形区域上的二重积分的数值计算
? 格式,
矩形区域的双重积分,
y=dblquad(Fun,xm,xM,ym,yM)
限定精度的双重积分,
y=dblquad(Fun,xm,xM,ym,yM,) ?
? 例,求解
>> f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');
>> y=dblquad(f,-2,2,-1,1)
y =
1.57449318974494
? 任意区域上二元函数的数值积分 (调用工具箱）
? 格式,
一般双重积分
J= quad2dggen(fun,ymin,ymax,xlower,xupper)
限定精度的双重积分
J= quad2dggen(fun,ymin,ymax,xlower,xupper,)
?
例
>> fh=inline('sqrt(1-x.^2/2)','x'); % 内积分上限
>> fl=inline('-sqrt(1-x.^2/2)','x'); % 内积分下限
>> f=inline('exp(-x.^2/2).*sin(x.^2+y)','y','x'); %
交换顺序的被积函数
>> y=quad2dggen(f,fl,fh,-1/2,1,eps)
y =
0.41192954617630
? 解析解方法,
>> syms x y
>> i1=int(exp(-x^2/2)*sin(x^2+y),y,-sqrt(1-x^2/2),sqrt(1-
x^2/2));
>> int(i1,x,-1/2,1)
Warning,Explicit integral could not be found,
> In D:\MATLAB6p5\toolbox\symbolic\@sym\int.m at line
58
ans =
int(2*exp(-1/2*x^2)*sin(x^2)*sin(1/2*(4-2*x^2)^(1/2)),x =
-1/2,,1)
>> vpa(ans)
ans =
.41192954617629511965175994017601
2 2
2
11
22
1 1
s i n ( )
y x
y
I e x y d x d y
?
?
? ??
????
2
22
22
1
sin ( )
x
xy
I e x y d x d y
?
??
????
例,计算单位圆域上的积分,
先把二重积分转化为二次积分的形式,
>> syms x y
i1=int(exp(-x^2/2)*sin(x^2+y),x,-sqrt(1-y.^2),sqrt(1-y.^2));
Warning,Explicit integral could not be found,
> In D:\MATLAB6p5\toolbox\symbolic\@sym\int.m at line 58
对 x是不可积的，故调用解析解方法不会得出结
果，而数值解求解不受此影响。
>> fh=inline('sqrt(1-y.^2)','y'); % 内积分上限
>> fl=inline('-sqrt(1-y.^2)','y'); % 内积分下限
>> f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y'); %交换顺序
的被积函数
>> I=quad2dggen(f,fl,fh,-1,1,eps)
Integral did not converge--singularity likely
I =
0.53686038269795
4.3.6 三重定积分的数值求解
? 格式,
I=triplequad(Fun,xm,xM,ym,yM,zm,zM，,@quadl)
其中 @quadl为具体求解一元积分的数值函数,
也可选用 @quad或自编积分函数,但格式要与
quadl一致。
?
? 例,
>> triplequad(inline('4*x.*z.*exp(-x.*x.*y-z.*z)',…
'x','y','z'),0,2,0,pi,0,pi,1e-7,@quadl)
ans =
3.10807945143834
