中科院的matlab课件,中科院matlab课件

最新推荐文章于 2022-05-06 10:00:00 发布

知行晓政

最新推荐文章于 2022-05-06 10:00:00 发布

阅读量443

点赞数

文章标签：中科院的matlab课件

PPT内容

这是中科院matlab课件，关于微积分问题的计算机求解，包括了微积分问题的解析解，函数的级数展开与级数求和问题求解，数值微分，数值积分问题，曲线积分与曲面积分的计算等内容，欢迎点击下载。

第三章微积分问题的计算机求解

微积分问题的解析解

函数的级数展开与级数求和问题求解

数值微分

数值积分问题

曲线积分与曲面积分的计算

3.1 微积分问题的解析解 3.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 =

在(-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)

3.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)

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)

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)

3.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)

------------

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)

例：

>> syms x; 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)

>> vpa(I1,70)

ans =

1.085853317666016569702419076542265042534236293532156326729917229308528

>> I2=int(exp(-x^2/2),x,0,inf)

I2 =

1/2*2^(1/2)*pi^(1/2)

多重积分问题的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) ???

f2=sin(x^2*y)/exp(y*x^2 + z^2)

>> simple(f1-f2)

ans =

顺序的改变使化简结果不同于原函数，但其误差为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*hypergeom([1],[2],-pi^2) ???

Ei(n,z)为指数积分，无解析解，但可求其数值解：

>> vpa(ans,60)

ans =

3.10807940208541272283461464767138521019142306317021863483588

3.2 函数的级数展开与级数求和问题求解

3.2.1 Taylor 幂级数展开

3.2.2 Fourier 级数展开

3.2.3 级数求和的计算

3.2.1 Taylor 幂级数展开 3.2.1.1 单变量函数的 Taylor 幂级数展开

例:

>> syms x; f=sin(x)/(x^2+4*x+3);

>> y1=taylor(f,x,9); pretty(y1)

23 34 4087 3067 515273 386459

1/3 x - 4/9 x2 + -- x3 - ---- x4 + ------x5 - ------ x6 +---------- x7 - --------- x8

54 81 9720 7290 1224720 918540

>> 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 +(-15697/6075000*sin(2)+28/50625*cos(2))*(x-2)^4 +(203/6075000*cos(2)+6277/11390625*sin(2))*(x-2)^5 +(-585671/2733750000*sin(2)-623/11390625*cos(2))*(x-2)^6 +(262453/19136250000*cos(2)+397361/5125781250*sin(2))*(x-2)^7 +(-875225059/34445250000000*sin(2)-131623/35880468750*cos(2))*(x-2)^8

>> syms a; 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)*a-4*sin(a)-2*sin(a)*a)/(a^2+3+4*a)^2*(4+2*a))/(a^2+3+4*a)*(x-a)^2+…

例:对y=sinx进行Taylor幂级数展开,并观察不同阶次的近似效果。

>> x0=-2*pi:0.01:2*pi; y0=sin(x0); syms x; 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 =

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/39916800*x^11

p =

x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/39916800*x^11+1/6227020800*x^13

p =

x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/39916800*x^11+1/6227020800*x^13-1/1307674368000*x^15

3.2.1.2 多变量函数的Taylor 幂级数展开

多变量函数在

的Taylor幂级数的展开

例：？？？

>> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y);

>> F=maple('mtaylor',f,'[x,y]',8)

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*x-y*x^3-y^2*x^2-2*y*x^4-3*y^2*x^3-2*y^3*x^2-y^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

>> syms a; 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*a-y)+(2*a-2)*exp(-a^2-y^2-a*y))*(x-a)+((a^2-2*a)*exp(-a^2-y^2-a*y)*(-1+2*a^2+2*a*y+1/2*y^2)+exp(-a^2-y^2-a*y)+(2*a-2)*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=-pi; b=pi; 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 ％换回变量区域

例:

>> syms x; f=x*(x-pi)*(x-2*pi);

>> [A,B,F]=fseries(f,x,6,0,2*pi)

A =

[ 0, 0, 0, 0, 0, 0, 0]

B =

[ -12, 3/2, -4/9, 3/16, -12/125, 1/18]

F =

12*sin(x)+3/2*sin(2*x)+4/9*sin(3*x)+3/16*sin(4*x)+12/125*sin(5*x)+1/18*sin(6*x)

例:

>> syms x; 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/pi*sin(x)+4/3/pi*sin(3*x)

……

3.2.3 级数求和的计算

是在符号工具箱中提供的

例:计算

>> format long; sum(2.^[0:63]) %数值计算

ans =

1.844674407370955e+019

>> sum(sym(2).^[0:200]) % 或 syms k; symsum(2^k,0,200)

％把2定义为符号量可使计算更精确

ans =

3213876088517980551083924184682325205044405987565585670602751

>> syms k; symsum(2^k,0,200)

ans =

3213876088517980551083924184682325205044405987565585670602751

例:试求解无穷级数的和

>> syms n; s=symsum(1/((3*n-2)*(3*n+1)),n,1,inf)

%采用符号运算工具箱

s =

1/3

>> m=1:10000000; s1=sum(1./((3*m-2).*(3*m+1)));%数值计算方法,双精度有效位16,“大数吃小数”,无法精确

>> format long; s1 % 以长型方式显示得出的结果

s1 =

0.33333332222165

例:求解

>> syms n x

>> s1=symsum(2/((2*n+1)*(2*x+1)^(2*n+1)),n, 0,inf);

>> simple(s1) % 对结果进行化简，MATLAB 6.5 及以前版本因本身 bug 化简很麻烦

ans =

log((((2*x+1)^2)^(1/2)+1)/(((2*x+1)^2)^(1/2)-1))

%实际应为log((x+1)/x)

例:求

>> syms m n; limit(symsum(1/m,m,1,n)-log(n),n,inf)

ans =

eulergamma

>> vpa(ans, 70) % 显示 70 位有效数字

ans =

.5772156649015328606065120900824024310421593359399235988057672348848677

符号函数计算器

单变量符号函数计算器

Taylor 逼近计算器

单变量符号函数计算器(1/3)

在命令窗口中执行 funtool 即可调出单变量符号函数计算器。单变量符号函数计算器用于对单变量函数进行操作，可以对符号函数进行化简、求导、绘制图形等。该工具的界面如图所示。

单变量符号函数计算器(2/3)

1．输入框的功能

如图：

单变量符号函数计算器(3/3)

单变量符号函数计算器应用示例

Taylor 逼近计算器

Taylor 逼近计算器用于实现函数的 taylor 逼近。在命令窗口中输入 taylortool，调出Taylor 逼近计算器，界面及功能如图。

MAPLE 函数的调用

maple 函数的使用

mfun 函数的使用

maple 函数的使用

maple 是符号工具箱中的一个通用命令，使用它可以实现对 MAPLE 中大部分函数的调用。其使用格式为：

r = maple('statement')，其中 statement 为符合 MAPLE 语法的可执行语句的字符串，该命令将 statement 传递给 MAPLE，该命令的输出结果也符合 MAPLE 的语法；

r = maple('function',arg1,arg2,...)，该函数调用引号中的函数，并接受指定的参数，相当于 MAPLE 语句 function(arg1,arg2,...)；

[r, status] = maple(...)，返回函数的运行状态，如果函数运行成功，则 status 为 0，r 为运行结果；如果函数运行失败，则 status 为一个正数，r 为相应的错误信息；

maple('traceon') 或者 maple trace on，输出 MAPLE 函数运行中的所有子表达式和运行结果；

maple('traceoff') 或 maple trace off，不显示中间过程。

mfun 函数的使用

mfun 函数用于对 maple 函数进行数字评估。该函数的调用格式为：

Y = mfun('function',par1,par2,par3,par4)。

该语句对指定的数学函数进行评估。其中 'function' 指定待评估的函数，par1、par2 等为 'function' 的参数，为待评估的数值，其维数有 'function' 函数的参数类型确定。在该语句中最多可以设置四个参数，最后一个参数可以为矩阵。

用户可以通过 help mfunlist 查看 MATLAB 中 mfun 可以调用的函数列表，另外，可以通过 mhelp function 查看指定函数的相关信息。

3.3.1 数值微分算法

向前差商公式：

向后差商公式

两种中心公式：

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;

调用格式：

y为等距实测数据， dy为得出的导数向量， dx为相应的自变量向量，dy、dx的数据比y短。

例：

求导数的解析解,再用数值微分求取原函数的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

3.3.3 用插值、拟合多项式的求导数

基本思想：当已知函数在一些离散点上的函数值时，该函数可用插值或拟合多项式来近似，然后对多项式进行微分求得导数。

选取x=0附近的少量点

进行多项式拟合或插值

g(x)在x=0处的k阶导数为

通过坐标变换用上述方法计算任意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=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

3.3.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/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])

3.4 数值积分问题 4.3.1 由给定数据进行梯形求积

格式： 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

>> format long; v

v =

0.10000000000000 0.05389175150076 0.01277491516591

0.01000000000000 0.06654169546584 0.00012497120083

0.00100000000000 0.06666541668004 0.00000124998663

0.00010000000000 0.06666665416667 0.00000001250000

0.00001000000000 0.06666666654167 0.00000000012500

0.00000100000000 0.06666666666542 0.00000000000125

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。后面函数算法更精确，精度更高。

例：

函数定义被积函数：

>> 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 =

.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

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

3.4.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)));

function y=gaussf(x)

y=cos(x);

>> gauss2('gaussf',0,1)

ans =

0.8415

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.57449318974494

任意区域上二元函数的数值积分 (调用工具箱NIT)，该函数指定顺序先x后y.

解析解方法：

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

例：计算单位圆域上的积分：

先把二重积分转化：

对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

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为曲线上某点的弧长,又称这类曲线积分为对弧长的曲线积分.

数学表示

若

弧长表示为

例：

>> syms t; syms a positive; 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)

％绘出两条曲线

>> syms x; 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))

3.5.1.2 第二类曲线积分

又称对坐标的曲线积分，起源于变力

沿曲线移动时作功的研究

曲线亦为向量,若曲线可以由参数方程表示

则两个向量的点乘可由这两个向量直接得出.

例：求曲线积分

>> syms t; syms a positive; 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

例：

>> syms x; 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 第一类曲面积分

其中为小区域的面积，故又称为对面积的曲面积分。曲面由给出，则该积分可转换成x-y平面的二重积分为

例：

％四个平面，其中三个被积函数的值为0，只须计算一个即可。

>> syms x y; syms a positive; 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

若曲面由参数方程

曲面积分

例：

>> syms u v; syms a positive;

>> 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].

>> syms u v; syms a b c positive;

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

中科院的matlab课件,中科院matlab课件

“相关推荐”对你有帮助么？