ImageVerifierCode 换一换
格式:PPT , 页数:69 ,大小:1.51MB ,
资源ID:501673      下载积分:10 文钱
快捷下载
登录下载
邮箱/手机:
温馨提示:
快捷下载时,用户名和密码都是您填写的邮箱或者手机号,方便查询和重复下载(系统自动生成)。 如填写123,账号就是123,密码也是123。
特别说明:
请自助下载,系统不会自动发送文件的哦; 如果您已付费,想二次下载,请登录后访问:我的下载记录
支付方式: 支付宝    微信支付   
验证码:   换一换

加入VIP,省得不是一点点
 

温馨提示:由于个人手机设置不同,如果发现不能下载,请复制以下地址【https://www.wenke99.com/d-501673.html】到电脑端继续下载(重复下载不扣费)。

已注册用户请登录:
账号:
密码:
验证码:   换一换
  忘记密码?
三方登录: QQ登录   微博登录 

下载须知

1: 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。
2: 试题试卷类文档,如果标题没有明确说明有答案则都视为没有答案,请知晓。
3: 文件的所有权益归上传用户所有。
4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
5. 本站仅提供交流平台,并不能对任何下载内容负责。
6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

版权提示 | 免责声明

本文(微积分问题的计算机求解.PPT)为本站会员(国***)主动上传,文客久久仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对上载内容本身不做任何修改或编辑。 若此文所含内容侵犯了您的版权或隐私,请立即通知文客久久(发送邮件至hr@wenke99.com或直接QQ联系客服),我们立即给予删除!

微积分问题的计算机求解.PPT

1、第四章 微积分问题的计算机求解,微积分问题的解析解函数的级数展开与级数求和问题求解*数值微分数值积分问题曲线积分与曲面积分的计算*,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(x3)-1)/(1-cos(sqrt(x-sin(x

2、),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不是确定的值,而是另一个

3、变量的函数,如x-g(y),则上述的极限求取顺序不能交换。,例:求出二元函数极限值 syms x y a; f=exp(-1/(y2+x2)*sin(x)2/x2*(1+1/y2)(x+a2*y2); L=limit(limit(f,x,1/sqrt(y),y,inf)L =exp(a2),4.1.2 函数导数的解析解,函数的导数和高阶导数格式: y=diff(fun,x) %求导数 y= diff(fun,x,n) %求n阶导数例: 一阶导数: syms x; f=sin(x)/(x2+4*x+3); f1=diff(f); pretty(f1),cos(x) sin(x) (2 x + 4

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); tocelapsed_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) (

5、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

6、(f,y,n),x,m)例: 求其偏导数并用图表示。 syms x y z=(x2-2*x)*exp(-x2-y2-x*y); zx=simple(diff(z,x)zx = -exp(-x2-y2-x*y)*(-2*x+2+2*x3+x2*y-4*x2-2*x*y), zy=diff(z,y)zy =(x2-2*x)*(-2*y-x)*exp(-x2-y2-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), con

7、tour(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(x2*y)*exp(-x2*y-z2); 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 ) (

8、cos(x y) - 10 cos(x y) y x + 4 2 4 2 2 4 2 2sin(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)

9、, 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=(x2-2*x)*exp(-x2-y2-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 - -

10、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)/(x2+4*x+3); y1=diff(y); y0=int(y1); pretty(y0) % 对导数积分 sin(x) sin(

11、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(x3*cos(a*x)2,x)f = 1/16*(4*a3*x3*sin(2*a*x)+2*a4 *x4+6*a2*x2*cos(2*a*x)-6*a*x*sin(2*a*x)-3*cos(2*a*x)-3)/a4 f1=x4/8+(x3/(4*a)-3*x/(8*a3)*sin

12、(2*a*x)+. (3*x2/(8*a2)-3/(16*a4)*cos(2*a*x); simple(f-f1) % 求两个结果的差ans = -3/16/a4,定积分与无穷积分计算:格式: I=int(f,x,a,b)格式: I=int(f,x,a,inf),多重积分问题的MATLAB求解例: syms x y z; f0=-4*z*exp(-x2*y-z2)*(cos(x2*y)-10*cos(x2*y)*y*x2+. 4*sin(x2*y)*x4*y2+4*cos(x2*y)*x4*y2-sin(x2*y); f1=int(f0,z);f1=int(f1,y);f1=int(f1,x)

13、; f1=simple(int(f1,x)f1 = exp(-x2*y-z2)*sin(x2*y), f2=int(f0,z); f2=int(f2,x); f2=int(f2,x); f2=simple(int(f2,y)f2 =2*exp(-x2*y-z2)*tan(1/2*x2*y)/(1+tan(1/2*x2*y)2) simple(f1-f2)ans =0 顺序的改变使化简结果不同于原函数,但其误差为0,表明二者实际完全一致。这是由于积分顺序不同,得不出实际的最简形式。,例: syms x y z int(int(int(4*x*z*exp(-x2*y-z2),x,0,1),y,0,

14、pi),z,0,pi)ans =(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2)*pi2*hypergeom(1,2,-pi2)Ei(n,z)为指数积分,无解析解,但可求其数值解: vpa(ans,60) ans = 3.10807940208541272283461464767138521019142306317021863483588,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

15、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*Dt2);L0=3;,case 3 dy=(-diff(yx1)+7*diff(yx2)-6*diff(yx3)-6*diff(yx4)+. 7*dif

16、f(yx5)-diff(yx6)/(8*Dt3); L0=5; case 4 dy = (-diff(yx1)+11*diff(yx2)-28*diff(yx3)+28* diff(yx4)-11*diff(yx5)+diff(yx6)/(6*Dt4); L0=5; end dy=dy(L0+1:end-L0); dx=(1:length(dy)+L0-2-(n2)*Dt;调用格式: y为 等距实测数据, dy为得出的导数向量, dx为相应的自变量向量 。,例:求导数的解析解,再用数值微分求取原函数的14 阶导数,并和解析解比较精度。 h=0.05; x=0:h:pi; syms x1; y=

17、sin(x1)/(x12+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)

18、,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阶导数为,通过坐标变换用上述

19、方法计算任意x点处的导数值令将g(x)写成z的表达式导数为可直接用 拟合节点 得到系数 d=polyfit(x-a,y,length(xd)-1),例:数据集合如下: 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=polyf

20、it(xd-a,yd,L-1);fact=1; for k=1:L-1;fact=factorial(k),fact;end deriv=d.*factderiv = 1.8750 -1.3750 1.0406 -0.9710 0.6533 0.6376,建立计算插值多项式各阶导数的poly_drv.mfunction 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);endder=c.*fact;例: xd= 0 0.200

21、0 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矩阵是建立在等间距的形式生成的网格基础上,则实际梯度为,例:计算梯度,绘制引力线图: 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

22、/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=

23、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),

24、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)/2S = 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; % 这样赋值能确保

25、 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 vv = 0.1000 0.0539 0.0128 0.0100 0.0665 0.0001 0.0010 0.0667 0.0000

26、 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, ) %限定精度的定积分求解,默认精度为106。,例:,第三种:匿名函数(MATLAB 7.0),第二种:inline 函数,第一种,一般函数方法,函数定义被积函数: y=quad(c3ffun,0,1.5)y = 0.9661用 inline 函数定义被积函数: f=inlin

27、e(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(-x2),0,1.5),60) y0 = .966105146475310713936933729949905794996224943257461473285749 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 forma

28、t 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).*(x2); y(end)=0; x=eps, x; y=0,y; fill(x,y,g),调用quad( ): f=inline(exp(x.2).*(x2)./(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;

29、 I=vpa(int(exp(x2),0,2)+int(80/(4-sin(16*pi*x),2,4) I = 57.764450125053010333315235385182,4.3.3 Gauss求积公式,为使求积公式得到较高的代数精度对求积区间a,b,通过变换有,以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.88888889*gaussf(x(2);funct

30、ion 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=(x2-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(s

31、p1,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=(x2-2*x)*exp(-x2-y2-x*y); ezsurf(diff(diff(z,x),y),-

32、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, 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

33、,-); 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) 限定精度的双重

34、积分 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(-x2/2)*sin(x2+y),y,-sqrt(1-x2/2),sqrt(1-x2/2)

35、; int(i1,x,-1/2,1)Warning: Explicit integral could not be found. In D:MATLAB6p5toolboxsymbolicsymint.m at line 58 ans = int(2*exp(-1/2*x2)*sin(x2)*sin(1/2*(4-2*x2)(1/2),x = -1/2 . 1) vpa(ans) ans = .41192954617629511965175994017601,例:计算单位圆域上的积分: 先把二重积分转化为二次积分的形式:, syms x y i1=int(exp(-x2/2)*sin(x2+y

36、),x,-sqrt(1-y.2),sqrt(1-y.2);Warning: Explicit integral could not be found. In D:MATLAB6p5toolboxsymbolicsymint.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

37、,-1,1,eps)Integral did not converge-singularity likelyI = 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,

Copyright © 2018-2021 Wenke99.com All rights reserved

工信部备案号浙ICP备20026746号-2  

公安局备案号:浙公网安备33038302330469号

本站为C2C交文档易平台,即用户上传的文档直接卖给下载用户,本站只是网络服务中间平台,所有原创文档下载所得归上传人所有,若您发现上传作品侵犯了您的权利,请立刻联系网站客服并提供证据,平台将在3个工作日内予以改正。