以下试题来自：
填空题由曲线L：绕y轴旋转一周所得到的旋转曲面在点处的指向外侧的单位法向量为_________． 正确答案：
为您推荐的考试题库
你可能感兴趣的试题
1.填空题 正确答案：2.填空题 正确答案：23.填空题 正确答案：4.填空题 正确答案：5.填空题 正确答案：7x一9y一5z+10=0
热门相关试卷
最新相关试卷matlab在科学计算中的应用3
来源：互联网
以下内容为系统自动转化的文字版，可能排版等有问题，仅供您参考：第三章 微积分问题的计算机求解? ? ? ? ? 微积分问题的解析解 函数的级数展开与级数求和问题求解 数值微分 数值积分问题 曲线积分与曲面积分的计算3.1 微积分问题的解析解3.1.1 极限问题的解析解? 单变量函数的极限– 格式1： L= limit( fun, x, x0)– 格式2： L= limit( fun, x, x0, ‘left’ 或 ‘right’)? 例: 试求解极限问题&& && f=x*(1+a/x)^x*sin(b/x); && L=limit(f,x,inf) L= exp(a)*b? 例:求解单边极限问题&& && 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)，则上述的极限求 取顺序不能交换。? 例：求出二元函数极限值&& && 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)3.1.2 函数导数的解析解? 函数的导数和高阶导数– 格式： y=diff(fun,x) %求导数 y= diff(fun,x,n) %求n阶导数? 例：一阶导数： && f=sin(x)/(x^2+4*x+3); && f1=diff(f); pretty(f1)cos(x) sin(x) (2 x + 4) --------------- - ------------------2 2 2 x +4x+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)? 例： 用图表示。求其偏导数并&& 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.*y4*x.^2-2*x.*y); && zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y); % 偏导 的数值解 && quiver(x,y,zx,zy) % 绘制引力线? 例&& 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 矩阵&& && x=r*sin(theta)*cos(phi); && y=r*sin(theta)*sin(phi); && z=r*cos(theta); && J=jacobian([x; 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)? 例：&& 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)3.1.3 积分问题的解析解? 不定积分的推导：–格式： F=int(fun,x)? 例： 用diff() 函数求其一阶导数，再积分，检验是否可以 得出一致的结果。&& 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? 例:证明&& 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)? 例： f ( x) ? e? x2 / 2&& I1=int(exp(-x^2/2),x,0,1.5) ％无解 I1 = 1/2*erf(3/4*2^(1/2))*2^(1/2)*pi^(1/2) 2 x ?t 2 erf ( x) ? e dt && vpa(I1,70) 0 ? ans = 1. && I2=int(exp(-x^2/2),x,0,inf) I2 = 1/2*2^(1/2)*pi^(1/2)?? 多重积分问题的MATLAB求解 ? 例：&& 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 %2009b && int(int(int(4*x*z*exp(-x^2*y-z^2),x,0,1),y,0,pi),z,0,pi) ans =-(1/exp(pi^2) - 1)*(eulergamma + log(pi) - Ei(-pi))Ei(n,z)为指数积分，无解析解，但可求其数值解： && vpa(ans,60) ans =1.3.2 函数的级数展开与 级数求和问题求解? 3.2.1 Taylor 幂级数展开? 3.2.2 Fourier 级数展开 ? 3.2.3 级数求和的计算3.2.1 Taylor 幂级数展开3.2.1.1 单变量函数的 Taylor 幂级数展开
例:&& f=sin(x)/(x^2+4*x+3); && y1=taylor(f,x,9); pretty(y1)2 23 3 34 4 7 6
1/3 x - 4/9 x + -- x - ---- x + ------x - ------ x +---------- x - --------- x 54 81
&& taylor(f,x,9,2)ans = 1/15*sin(2)+(1/15*cos(2)-8/225*sin(2))*(x-2)+ (-127/6750*sin(2)8/225*cos(2))*(x-2)^2 +(23/6750*cos(2)+628/50625*sin(2))*(x-2)^3 +(0*sin(2)+28/50625*cos(2))*(x-2)^4 +(203/6075000*cos(2)+5*sin(2))*(x-2)^5 +(3750000*sin(2)-623/*cos(2))*(x-2)^6 +(*cos(2)+5781250*sin(2))*(x-2)^7 +(*sin(2)-*cos(2))*(x-2)^8&& taylor(f,x,5,a) % 结果较冗长，显示从略ans = sin(a)/(a^2+3+4*a) +(cos(a)sin(a)/(a^2+3+4*a)*(4+2*a))/(a^2+3+4*a)*(x-a) +(sin(a)/(a^2+3+4*a)-1/2*sin(a)-(cos(a)*a^2+3*cos(a)+4*cos(a)*a4*sin(a)-2*sin(a)*a)/(a^2+3+4*a)^2*(4+2*a))/(a^2+3+4*a)*(xa)^2+…例:对y=sinx进行Taylor幂级数展开,并观察不同阶次的近 似效果。 && x0=-2*pi:0.01:2* y0=sin(x0); y=sin(x); && plot(x0,y0,'r-.'), axis([-2*pi,2*pi,-1.5,1.5]); hold on && for n=[8:2:16] p=taylor(y,x,n), y1=subs(p,x,x0); line(x0,y1) end % p=taylor(y,n) (2009b)p= x-1/6*x^3+1/120*x^5-1/5040*x^7 p= x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9 p= x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/*x^11 p= x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^91/*x^11+1/*x^13p= x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^91/*x^11+1/*x^131/0*x^153.2.1.2 多变量函数的Taylor 幂级数展开? 多变量函数 f ( x1 , x2 ,?, xn ) 在 (a1, a2 ,?, an ) 的Taylor幂级数的展开? 例：？？？&& f=(x^2-2*x)*exp(-x^2-y^2-x*y); && F=maple('mtaylor',f,'[x,y]',8) % Maple engine is notinstalled.(2009b)F= mtaylor((x^2-2*x)*exp(-x^2-y^2-x*y),[x, y],8)&& maple(‘readlib(mtaylor)’);％读库，把函数调入内存 && F=maple('mtaylor',f,'[x,y]',8)F= -2*x+x^2+2*x^3-x^4-x^5+1/2*x^6+1/3*x^7+2*y*x^2+2*y^2*xy*x^3-y^2*x^2-2*y*x^4-3*y^2*x^3-2*y^3*x^2y^4*x+y*x^5+3/2*y^2*x^4+y^3*x^3+1/2*y^4*x^2+y*x^6+2*y^ 2*x^5+7/3*y^3*x^4+2*y^4*x^3+y^5*x^2+1/3*y^6*x&& F=maple('mtaylor',f,'[x=1,y=a]',3); && F=maple('mtaylor',f,'[x=a]',3)F= (a^2-2*a)*exp(-a^2-y^2-a*y)+((a^2-2*a)*exp(-a^2-y^2-a*y)*(-2*ay)+(2*a-2)*exp(-a^2-y^2-a*y))*(x-a)+((a^2-2*a)*exp(-a^2-y^2a*y)*(-1+2*a^2+2*a*y+1/2*y^2)+exp(-a^2-y^2-a*y)+(2*a2)*exp(-a^2-y^2-a*y)*(-2*a-y))*(x-a)^23.2.2 Fourier 级数展开function [A,B,F]=fseries(f,x,n,a,b) if nargin==3, a=- b= end L=(b-a)/2; if a+b, f=subs(f,x,x-L-a); end％变量区域互换 A=int(f,x,-L,L)/L; B=[]; F=A/2; %计算a0 for i=1:n an=int(f*cos(i*pi*x/L),x,-L,L)/L; bn=int(f*sin(i*pi*x/L),x,-L,L)/L; A=[A, an]; B=[B,bn]; F=F+an*cos(i*pi*x/L)+bn*sin(i*pi*x/L); end if a+b, F=subs(F,x,x+L+a); end ％换回变量区域例:&& f=x*(x-pi)*(x-2*pi); && [A,B,F]=fseries(f,x,6,0,2*pi)A= [ (-16)*pi^3, 24*pi, (-6)*pi, (8*pi)/3, -(3*pi)/2, (24*pi)/25, -(2*pi)/3] B= [ -(12*pi - 24*pi^3)/pi, ((3*pi)/2 - 12*pi^3)/pi, -((4*pi)/9 - 8*pi^3)/pi, ((3*pi)/16 6*pi^3)/pi, (24*pi^2)/5 - 12/125, (pi/18 - 4*pi^3)/pi] F= (sin(x)*(12*pi - 24*pi^3))/pi - 24*pi*cos(x) - 8*pi^3 - 6*pi*cos(2*x) (8*pi*cos(3*x))/3 - (3*pi*cos(4*x))/2 - (24*pi*cos(5*x))/25 (2*pi*cos(6*x))/3 - sin(5*x)*((24*pi^2)/5 - 12/125) + (sin(2*x)*((3*pi)/2 12*pi^3))/pi + (sin(3*x)*((4*pi)/9 - 8*pi^3))/pi + (sin(4*x)*((3*pi)/16 6*pi^3))/pi + (sin(6*x)*(pi/18 - 4*pi^3))/pi例:&& f=abs(x)/x; % 定义方波信号 && xx=[-pi:pi/200:pi]; xx=xx(xx~=0); xx=sort([xx,eps,eps]); % 剔除零点 && yy=subs(f,x,xx); plot(xx,yy,'r-.'), hold on % 绘 制出理论值并保持坐标系 && for n=2:20 [a,b,f1]=fseries(f,x,n), y1=subs(f1,x,xx); plot(xx,y1) enda= [ 0, 0, 0] b= [ 4/pi, 0] f1 = 4/pi*sin(x) a= [ 0, 0, 0, 0] b= [ 4/pi, 0, 4/(3*pi)] f1 = (4*sin(3*x))/(3*pi) + (4*sin(x))/pi……3.2.3 级数求和的计算? 是在符号工具箱中提供的例:计算&& sum(2.^[0:63]) %数值计算ans = 1.955e+019&& sum(sym(2).^[0:200]) % 或 symsum(2^k,0,200) ％把2定义为符号量可使计算更精确ans = 51&& symsum(2^k,0,200)ans = 51例:试求解无穷级数的和&& s=symsum(1/((3*n-2)*(3*n+1)),n,1,inf) %采用符号运算工具箱 s= 1/3 && m=1:; s1=sum(1./((3*m-2).*(3*m+1)));%数 值计算方法,双精度有效位16,“大数吃小数”,无法精 确 && s1 % 以长型方式显示得出的结果 s1 = 0.65例:求解 && syms n x && s1=symsum(2/((2*n+1)*(2*x+1)^(2*n+1)),n, 0,inf); && simple(s1) % 对结果进行化简，MATLAB 6.5 及以前版本因本身 bug 化简很麻烦. 2009b??? ans = log((((2*x+1)^2)^(1/2)+1)/(((2*x+1)^2)^(1/2)-1)) %实际应为log((x+1)/x)例:求&& limit(symsum(1/m,m,1,n)-log(n),n,inf) ans = eulergamma && vpa(ans, 70) % 显示 70 位有效数字 ans = .73.3 数值微分差商型求导公式 f ( x ? h) ? f ( x ) 由导数定义 f ( x) ? lim h ?0 h'（ 1 ）向前差商公式'f ( x ? h) ? f ( x ) f ( x) ? h （2）向后差商公式 f ( x ) ? f ( x ? h) ' f ( x) ? h （3）中心差商公式 （中点方法 ） f ( x ? h) ? f ( x ? h) f ( x) ? 2h'A f(x) CTBx-hxx+h3.3.1 数值微分算法? 向前差商公式： ? 向后差商公式两种中心公式：' 2 '' 3 ''' 4 f ( x ) ? ? tf ( x ) ? ? t f ( x ) / 2! ? ? t f ( ? ) / 3! ? o ( ? t ) % f ( x) ? 2?t f ( x) ? ?tf ' ( x) ? ?t 2 f '' ( x) / 2!? ?t 3 f ''' (? ) / 3!? o( ?t 4 ) ? 2?t 2 ? t ? f ' ( x) ? f ''' (? ) 3!
3.3.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; ％数值计算diff(X)表示数组X相邻两数的差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; 调用格式：[d y , dx ] ? diff _ ctr ( y, ?t, n)y为 等距实测数据， dy为得出的导数向量， dx为 相应的自变量向量，dy、dx的数据比y短 。? 例： 求导数的解析解,再用数值微分求取原函数的1～ 4 阶导数，并和解析解比较精度。&& h=0.05; x=0:h: && 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.3.3.3 用插值、拟合多项式的求导数? 基本思想：当已知函数在一些离散点上的函 数值时，该函数可用插值或拟合多项式来近 似，然后对多项式进行微分求得导数。 ? 选取x=0附近的少量点( xi , yi ), i ? 1, 2,?, n ? 1? 进行多项式拟合或插值g( x) ? c1x ? c2 x ??? cn x ? cn?1nn?1? g(x)在x=0处的k阶导数为g ( k ) (0) ? cn?1?k k ! k ? 0,1, 2,?, n? 通过坐标变换用上述方法计算任意x点处的导 数值 z ? x?a ? 令 ? 将g(x)写成z的表达式? 导数为g ( x) ? g ( z) ? d1z ? d2 z ??? dn z ? dn?1nn?1g (a) ? g (0) ? dn?1?k k ! k ? 0,1,?, n(k ) (k )? 可直接用 g ( z ) 拟合节点 ( xi ? a, yi ) 得到系数 d i d=polyfit(x-a,y,length(xd)-1) d=polyfit(x-a,y,length(x)-1)? 例：数据集合如下： xd: 0 0.0 0.0 1.000 yd: 0.2 0.1 0.3 计算x=a=0.3处的各阶导数。 && xd=[ 0 0.0 0.0 1.000]; && yd=[0.2 0.1 0.3]; && 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.0 1.0 0.6? 建立用拟合（插值）多项式计算各阶导数的 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.* ? 例： && xd=[ 0 0.0 0.0 1.000]; && yd=[0.2 0.1 0.3]; && a=0.3; der=poly_drv(xd,yd,a) der = 0.0 1.0 1.87503.3.4 二元函数的梯度计算z ? f ( x, y)? 格式：[ f x , f y ] ? gradient ( z)? 若z矩阵是建立在等间距的形式生成的网格基 础上，则实际梯度为f x ? f x / ?x,f y ? f y / ?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); && && quiver(x,y,fx,fy) %绘制等高线与 引力线图? 绘制误差曲面:&& zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-4*x.^22*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]) && 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]) && surf(x,y,abs(fy-zy)); axis([-3 3 -2 2 0,0.06])3.4 数值积分问题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.0 1.9995 && S1=trapz(x1,y) % 得出和上述完全一致的结果 S1 = 1.0 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.01,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.00 0.76 0.910.00 0.00 0.00 0.00 0.000.84 0.04 0.67 0.67 0.420.83 0.63 0.00 0.00 0.253.4.2 单变量数值积分问题求解梯形公式? 格式：(变步长）(Fun:函数的字符串变量） y=quad(Fun,a,b) y=quadl(Fun,a,b) % 求定积分 y=quad(Fun,a,b, ? ) y=quadl(Fun,a,b,? ) %限定精 度的定积分求解，默认精度为10－6。后面函数 算法更精确，精度更高。? 例：第一种，一般函数方法第二种：inline 函数第三种：匿名函数(MATLAB 7.0)函数定义被积函数： && 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 = .9 && y=quad(f,0,1.5,1e-20) % 设置高精度，但该方法失效提高求解精度： && y=quadl(f,0,1.5,1e-20) y = 0.9661 && abs(y-y0) ans = .3892e-16 && format long ％16位精度 && y=quadl(f,0,1.5,1e-20) y = 0.31? 例：求解绘制函数:&& 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.63? 调用quadl( ):&& I2=quadl(f,0,4) I2 = 57.68&& I=vpa(int(exp(x^2),0,2)+int(80/(4-sin(16*pi*x)),2,4)) I= 57.3851823.4.3 Gauss求积公式? 为使求积公式得到较高的代数精度?1?1f ( x)dx ? ? Ak f ( xk )k ?0nn ? 1, x0 ? x1 ? ?0.5773503, A0 ? A1 ? 1.; n ? 2, x0 ? x2 ? ?0.7745967, A0 ? A2 ? 0. ? 0., A1 ? 0.;? 对求积区间[a,b]，通过变换 x ? b ? a t ? b ? a 2 2 ? 有?bab?a 1 b?a a ?b b?a n b?a a ?b f ( x)dx ? f( t? )dt ? Ak f ( t? ) ? ? 2 ?1 2 2 2 k ?0 2 2b?a b?a ? 以n=2的高斯公式为例： x? t? 2 2 function g=gauss2(fun,a,b) n ? 2, x0 ? x2 ? ?0.7745967, x1 ? 0. h=(b-a)/2; A0 ? A2 ? 0., A1 ? 0.; c=(a+b)/2; x=[h*(-0.7745967)+c, c, h*0.7745967+c]; g=h*(0.*(gaussf(x(1))+gaussf(x(3)))+0.* gaussf(x(2)));function y=gaussf(x) y=cos(x); && gauss2('gaussf',0,1) ans = 0.8415?bab?a n b?a a ?b f ( x)dx ? A f ( t? ) ? k 2 k ?0 2 23.4.4 双重积分问题的数值解? 矩形区域上的二重积分的数值计算? 格式： 矩形区域的双重积分： 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.94? 任意区域上二元函数的数值积分 (调用工具箱 NIT），该函数指定顺序先x后y.例&& 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.30? 解析解方法： && syms x y && i1=int(exp(-x^2/2)*sin(x^2+y), y, -sqrt(1-x^2/2), sqrt(1x^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 = .例：计算单位圆域上的积分：I ?x 2 ? y 2 ?1??ex2 ? 2sin( x 2 ? y )dxdy先把二重积分转化：I?&& syms x yi1=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?1?11? y 2? 1? y 2?ex2 ? 2sin( x 2 ? y )dxdy对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.953.4.5 三重定积分的数值求解? 格式: 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, 1, 0, pi, 0, pi,1e-7,@quadl)ans = 1.73283.5 曲线积分与曲面积分的计算? 3.5.1 曲线积分及MATLAB求解–第一类曲线积分 起源于对不均匀分布的空间曲线总 质量的求取.设空间曲线L的密度函数为 f(x,y,z),则其总质量其中s为曲线上某点的弧长,又称这类 曲线积分为对弧长的曲线积分.? 数学表示 若 x ? x(t ), y ? y(t ), z ? z (t )弧长表示为例：&& x=a*cos(t); y=a*sin(t); z=a*t;&& I=int(z^2/(x^2+y^2)*sqrt(diff(x,t)^2+diff(y,t)^2+ diff(z,t)^2),t,0,2*pi) I= 8/3*pi^3*a*2^(1/2)&& pretty(I)3 1/2 8/3 pi a 2例：&& x=0:.001:1.2; y1=x; y2=x.^2; plot(x,y1,x,y2) ％绘出两条曲线 && y1=x; y2=x^2; I1=int((x^2+y2^2)*sqrt(1+diff(y2,x)^2),x,0,1); && I2=int((x^2+y1^2)*sqrt(1+diff(y1,x)^2),x,1,0); I=I2+I1 I= -2/3*2^(1/2)+349/768*5^(1/2)+7/512*log(-2+5^(1/2))r 又称对坐标的曲线积分，起源于变力 f ( x, y, z) 沿曲线 l 移动时作功的研究3.5.1.2 第二类曲线积分r 曲线 ds 亦为向量,若曲线可以由参数方程表示则两个向量的点乘可由这两个向量直接得出.例：求曲线积分?lx? y x? y dx ? 2 dy 2 2 2 x ?y x ?y&& x=a*cos(t); y=a*sin(t); && F=[(x+y)/(x^2+y^2),-(x-y)/(x^2+y^2)]; ds=[diff(x,t);diff(y,t)]; && I=int(F*ds,t,2*pi,0) % 正向圆周 I= 2*pi例：&& y=x^2; F=[x^2-2*x*y,y^2-2*x*y]; ds=[1; diff(y,x)]; && I=int(F*ds,x,-1,1) I= -14/153.5.2曲面积分与MATLAB语言求解3.5.2.1 第一类曲面积分其中 dS 为小区域的面积，故又称为对面积 的曲面积分。曲面 S 由 z ? f ( x, y) 给出，则该积 分可转换成x-y平面的二重积分为例：％四个平面，其中三个被积函数的值为0，只 须计算一个即可。 && z=a-x-y; && I=int(int(x*y*z*sqrt(1+diff(z,x)^2+ diff(z,y)^2),y,0,a-x),x,0,a)I= 1/120*3^(1/2)*a^5若曲面由参数方程曲面积分例：&& && x=u*cos(v); y=u*sin(v); z=v;f=x^2*y+z*y^2; && E=simple(diff(x,u)^2+diff(y,u)^2+diff(z,u)^2); && F=diff(x,u)*diff(x,v)+diff(y,u)*diff(y,v)+diff(z,u)* diff(z,v); && G=simple(diff(x,v)^2+diff(y,v)^2+diff(z,v)^2); && I=int(int(f*sqrt(E*G-F^2),u,0,a),v,0,2*pi) I= 1/4*a*(a^2+1)^(3/2)*pi^2+1/8*log(-a+(a^2+1)^(1/2)) *pi^2-1/8*(a^2+1)^(1/2)*a*pi^23.5.2.2 第二类曲面积分又称对坐标的曲面积分可转化成第一类曲面积分
若曲面由参数方程给出例：的上半部，且积分沿椭球面的上面。 ％引入参数方程 x=a*sin(u)*cos(v); y=b*sin(u)*sin(v); z=c*cos(u)， u[0,pi/2], v[0,2*pi]. && && x=a*sin(u)*cos(v); y=b*sin(u)*sin(v); z=c*cos(u); && A=diff(y,u)*diff(z,v)-diff(z,u)*diff(y,v); && I=int(int(x^3*A,u,0,pi/2),v,0,2*pi) I = 2/5*pi*a^3*c*b[]相关文章： 1/21页