函数的数值导数和切平面
数学实验第五次课 微积分实验
4.1 函数的数值导数和切平面
4.1.1 法线
【例4.5.1-1】曲面法线演示。
y=-1:0.1:1;x=2*cos(asin(y)); [X,Y,Z]=cylinder(x,20);
surfnorm(X(:,11:21),Y(:,11:21),Z(:,11:21)); view([120,18])
图 4.5.1-1
4.1.2 偏导数和梯度 4.1.2.1 理论定义 4.1.2.2 数值计算指令
【例4.5.2.2-1】用一个简单矩阵表现diff 和gradient 指令计算方式。
F=[1,2,3;4,5,6;7,8,9] Dx=diff(F) Dx_2=diff(F,1,2) [FX,FY]=gradient(F)
[FX_2,FY_2]=gradient(F,0.5)
【例 4.5.2.2-2】研究偶极子(Dipole)的电势(Electric potential)和电场强度(Electric field density)。设在(a , b )
处有电荷+q ,在(-a , -b ) 处有电荷-q 。那么在电荷所在平面上任何一点的电势和场强分别为
q 11
,E =-∇V 。其中22
22
。
V (x , y ) =4πε
(
r +=
(x -a ) +(y -b ) , r -=
(x +a ) +(y +b )
r -
+
r ) -
1=9⋅109。又设电荷q =2⋅10-6
4πε
,a =1. 5,b =-1. 5。
clear;clf;q=2e-6;k=9e9;a=1.5;b=-1.5;x=-6:0.6:6;y=x; [X,Y]=meshgrid(x,y);
rp=sqrt((X-a).^2+(Y-b).^2);rm=sqrt((X+a).^2+(Y+b).^2); V=q*k*(1./rp-1./rm); [Ex,Ey]=gradient(-V);
AE=sqrt(Ex.^2+Ey.^2);Ex=Ex./AE;Ey=Ey./AE; cv=linspace(min(min(V)),max(max(V)),49); contourf(X,Y,V,cv,'k-') %axis('square')
title('\fontname{隶书}\fontsize{22}偶极子的场'),hold on quiver(X,Y,Ex,Ey,0.7) plot(a,b,'wo',a,b,'w+') plot(-a,-b,'wo',-a,-b,'w-')
xlabel('x');ylabel('y'),hold off
图 4.5.2.2-1
4.2 函数的零点
4.2.1 多项式的根 4.2.2 一元函数的零点
4.2.2.1 利用MATLAB 作图指令获取初步近似解
4.2.2.2 任意一元函数零点的精确解
【例 4.6.2.2-1】通过求f (t ) =(sin(1)
y=inline('sin(t)^2*exp(-a*t)-b*abs(t)','t','a','b');
%
2
t ) e
-at
-b t 的零点,综合叙述相关指令的用法。
(2)
a=0.1;b=0.5;t=-10:0.01:10; y_char=vectorize(y); Y=feval(y_char,t,a,b);
clf,plot(t,Y,'r');hold on,plot(t,zeros(size(t)),'k'); xlabel('t');ylabel('y(t)'),hold off
%
图4.6-1
(3) 由于Notebook 中无法实现zoom 、ginput 指令涉及的图形和鼠标交互操作,因此下面指令必须在MA TLAB 指令窗中运行,并得到如图4.6-2所示的局部放大图及鼠标操作线。
zoom on
[tt,yy]=ginput(5);zoom off
图 4.6-2
tt
(4)
[t4,y4,exitflag]=fzero(y,tt(4),[],a,b) %
(5)
[t3,y3,exitflag]=fzero(y,tt(3),[],a,b)
(6)
op=optimset('fzero')
op=optimset('tolx',0.01); op.TolX ans = 0.0100
(7)
[t4n,y4n,exitflag]=fzero(y,tt(4),op,a,b)
4.2.3 多元函数的零点
【例 4.6.3-1】求解二元函数方程组⎧f ⎨1(x , y ) =sin(x -y ) =0
⎩f 2
(x , y ) =cos(x +y ) =0的零点。
(1)
x=-2:0.5:2;y=x;[X,Y]=meshgrid(x,y); F1=sin(X-Y);F2=cos(X+Y); v=[-0.2, 0, 0.2]; contour(X,Y,F1,v)
hold on,contour(X,Y,F2,v),hold off
图4.6-3
(2)
[x0,y0]=ginput(2); disp([x0,y0]) -0.7926 -0.7843
0.7926 0.7843
(3)
fun='[sin(x(1)-x(2)),cos(x(1)+x(2))]'; [xy,f,exit]=fsolve(fun,[x0(2),y0(2)]) Optimization terminated successfully:
First-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected xy =
0.7854 0.7854 f =
1.0e-006 *
-0.0984 0.2011 exit = 1
% %
〖说明〗
[fun.m]
function ff=fun(x) ff(1)=sin(x(1)-x(2)); ff(2)=cos(x(1)+x(2));
4.3 函数极值点
4.3.1 一元函数的极小值点 4.3.2 多元函数的极小值点
【例4.7.2-1】求f (x , y ) =100(y -x ) +(1-x ) 的极小值点。它即是著名的Rosenbrock's "Banana" 测试函数。该测试函数有一片浅谷,许多算法难以越过此谷。(演示本例搜索过程的文件名为exm04072_1_1.m 。) (1)
ff=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x');
2
2
2
(2)
x0=[-1.2,1];[sx,sfval,sexit,soutput]=fminsearch(ff,x0) sx =
1.0000 1.0000 sfval = 8.1777e-010 sexit = 1 soutput =
iterations: 85 funcCount: 159
algorithm: 'Nelder-Mead simplex direct search'
(3)
[ux,sfval,uexit,uoutput,grid,hess]=fminunc(ff,x0)
Warning: Gradient must be provided for trust-region method; using line-search method instead.
> In D:\MATLAB6P1\toolbox\optim\fminunc.m at line 211
Optimization terminated successfully:
Current search direction is a descent direction, and magnitude of directional derivative in search direction less than 2*options.TolFun ux =
1.0000 1.0000 sfval = 1.9116e-011 uexit = 1 uoutput =
iterations: 26 funcCount: 162 stepsize: 1.2992 firstorderopt: 5.0020e-004
algorithm: 'medium-scale: Quasi-Newton line search' grid = 1.0e-003 * -0.5002 -0.1888 hess =
820.4028 -409.5496 -409.5496 204.7720
4.4 数值积分
4.4.1 一元函数的数值积分 4.4.1.1 闭型数值积分
【例 4.8.1.1-1】求I =⎰
1
e
-x
2
dx ,其精确值为0. 74684204 。
(1)
syms x;IS=int('exp(-x*x)','x',0,1) vpa(IS) IS =
1/2*erf(1)*pi^(1/2) ans =
.[***********][1**********]185
(2)
fun=inline('exp(-x.*x)','x');
Isim=quad(fun,0,1),IL=quadl(fun,0,1) Isim = 0.7468 IL = 0.7468
(3)
Ig=gauss10(fun,0,1) Ig =
0.7463
(4)
xx=0:0.1:1.5;ff=exp(-xx.^2); pp=spline(xx,ff); int_pp=fnint(pp);
Ssp=ppval(int_pp,[0,1])*[-1;1] Ssp = 0.7468
(5)
图4.8-1
4.4.1.2 开型数值积分
[gauss10.m]
function g = gauss10(fun,a,b) %GAUSS10(fun,a,b)
% fun
%====================================================== x = [0.1488743390;0.4333953941;0.6974095683;... 0.8650633667;0.9739065285];
w = [0.2955242247;0.2692667193;0.2190863625;... 0.1494513492;0.0666713443]; t = .5*(b+a)+.5*(b-a)*[-flipud(x);x]; W = [flipud(w);w];
g = sum(W.*feval(fun,t))*(b-a)/2;
【例 4.8.1.2-1】当f (x ) =cos(x ) 时,比较解析积分和近似积分。 (1)
syms x;F=int('cos(x)','x',-1,1),vpa(F) F = 2*sin(1) ans =
1.[***********][1**********]06 (2)
aF=cos(1/sqrt(3))+cos(-1/sqrt(3)) aF =
1.6758
【例4.8.1.2-2】求I =⎰
1
10
ln
x
dx ,准确结果是
2
=.88622692 。
(1)
syms x;IS=vpa(int('sqrt(log(1/x))','x',0,1)) Warning: Explicit integral could not be found.
> In D:\MATLAB6P5\toolbox\symbolic\@sym\int.m at line 58 In D:\MATLAB6P5\toolbox\symbolic\@char\int.m at line 9 IS =
.[***********][1**********]057
(2)用quad 指令求积
ff=inline('sqrt(log(1./x))','x');Isim=quad(ff,0,1) Warning: Divide by zero.
> In D:\MATLAB6P5\toolbox\matlab\funfun\inlineeval.m at line 13 In D:\MATLAB6P5\toolbox\matlab\funfun\@inline\feval.m at line 34 In D:\MATLAB6P5\toolbox\matlab\funfun\quad.m at line 62 Isim = 0.8862
(3)
Ig=gauss10(ff,0,1) Ig =
0.8861
4.4.2 多重数值积分
4.4.2.1 积分限为常数的二重积分指令
【例 4.8.2.1-1】计算S x 01=⎰
2
⎡11
x y dx ⎤dy 和S ⎰1
⎡2⎢⎣⎰0⎥⎦x 12=0⎢⎣⎰1x y dx ⎤⎥⎦
dy 。(1)
syms x y
ssx01=vpa(int(int(x^y,x,0,1),y,1,2)) ssx12=vpa(int(int(x^y,y,0,1),x,1,2))
Warning: Explicit integral could not be found.
> In D:\MATLAB6P5\toolbox\symbolic\@sym\int.m at line 58 ssx01 =
.[***********][1**********]435 ssx12 =
1.[***********][1**********]15
(2)
zz=inline('x.^y','x','y'); nsx01=dblquad(zz,0,1,1,2) nsx12=dblquad(zz,1,2,0,1) nsx01 = 0.4055 nsx12 =
1.2293
4.4.2.2 内积分限为函数的二重积分
[double_int.m]
function SS=double_int(fun,innlow,innhi,outlow,outhi) %double_int %
%fun %innlow,innhi
%outlow,outhi
y1=outlow;y2=outhi;x1=innlow;x2=innhi;f_p=fun; SS=quad(@G_yi,y1,y2,[],[],x1,x2,f_p);
[G_yi.m]
function f=G_yi(y,x1,x2,f_p) %G_yi %y %x1,x2 %
%f_p
y=y(:);n=length(y);
if isnumeric(x1)==1;xx1=x1*ones(size(y));else xx1=feval(x1,y);end if isnumeric(x2)==1;xx2=x2*ones(size(y));else xx2=feval(x2,y);end for i=1:n
f(i)=quad(f_p,xx1(i),xx2(i),[],[],y(i)); end
f=f(:);
【例 4.8.2.2-1】计算I =⎰
4
⎡21
⎢⎣(x 2+y 2y ) dx ⎤⎥⎦
dy 。 (1)
[x_low.m]
function f=x_low(y) f=sqrt(y);
(2)
(3)
ff=inline('x.^2+y.^2','x','y'); SS=double_int(ff,@x_low,2,1,4)
Warning: Minimum step size reached; singularity possible. > In D:\MATLAB6P5\toolbox\matlab\funfun\quad.m at line 88 In D:\MATLAB6p5\work\G_yi.m at line 11
In D:\MATLAB6P5\toolbox\matlab\funfun\@inline\feval.m at line 20 In D:\MATLAB6P5\toolbox\matlab\funfun\quad.m at line 62 In D:\MATLAB6p5\work\double_int.m at line 8 SS =
9.5810
(4)
Ssym=vpa(int(int('x^2+y^2','x','sqrt(y)',2),'y',1,4)) Ssym =
9.[***********][1**********]24