MATLAB习题答案(清华大学)_第1页
MATLAB习题答案(清华大学)_第2页
MATLAB习题答案(清华大学)_第3页
MATLAB习题答案(清华大学)_第4页
MATLAB习题答案(清华大学)_第5页
已阅读5页，还剩281页未读， 继续免费阅读

下载本文档

版权说明：本文档由用户提供并上传，收益归属内容提供方，若内容存在侵权，请进行举报或认领

文档简介

MATLAB习题答案（清华大学）

高等应用数学问题MATLAB求解

习题参考解答（薛定宇著）

目录

第1章计算机数学语言概述2

第2章MATLAB语言程序设计基础5

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

第4章线性代数问题的计算机求解29

第5章积分变换与复变函数问题的计算机求解43

第6章代数方程与最优化问题的计算机求解53

第7章微分方程问题的计算机求解71

第8章数据插值、函数逼近问题的计算机求解93

第9章概率论与数理统计问题的计算机求解114

第10章数学问题的非传统解法127

第A章自由数学语言Scilab简介136

第1章计算机数学语言概述

1在你的机器上安装MATLAB语言环境，并键入demo命令，由给出的菜单系统和对话

框原型

演示程序，领略MATLAB语言在求解数学问题方面的能力与方法。

【求解】在MATLAB提示符>>下键入demo命令，则将打开如图1-1所示的窗口，窗口

左侧

的列表框可以选择各种不同组合的演示内容。

图ITMATLAB演示程序界面1

例如，用户选择MATLAB!Graphics!VolumeVlsulization演示，则将得出如图1-2

所示的

演示说明，单击其中的Runthisdemo栏目，则将得出如图1-3所示的演示界面。用

户可以在

该界面下按按钮，逐步演示相关内容，而实现这样演示的语句将在该程序界面的下部窗

口中

给出。

2作者用MATLAB语言编写了给出例子的源程序，读者可以自己用type语句阅读一下

源程

序，对照数学问题初步理解语句的含义，编写的源程序说明由下表列出。

第1章计算机数学语言概述3

图1-2MATLAB演示程序界面举例

序号文件名程序说明

例1.1clexl.m利用MATLAB的符号运算工具箱求解微分问题

例L2clex2.m分别利用MATLAB的符号运算工具箱和数值运算功能求解多项式方程，

其中用数值方法得出

的结果有误差

例1.3clex3.m分别利用MATLAB的符号运算工具箱和数值运算功能计算Hilbert矩阵

的行列式，其中用数值

方法得出的结果有很大误差

例1.4clex4.m令xl=y;x2=y_,则可以将原来的二阶微分方程转换成2一阶微分

方程组，然后就可以求解微分

方程的数值解了，原方程是非线性微分方程，故不存在解析解。。de45()函数可以求解

常微分方

程组，而dde23()可以求解延迟微分方程，或更直观地采用Simulink绘制求解框图。

例1.5clex5.m线性规划问题调用最优化工具箱中的1inprog0函数可以立即得出结

果，若想求解整数规划问

题，则需要首先安装整数规划程序ipslvmex(),

4第1章计算机数学语言概述

图1-3MATLAB体视化演示程序界面

第2章MATLAB语言程序设计基础

1启动MATLAB环境,并给出语句tic,Agrand(500);B=inv(A);norm(A*B-eye(500)),

toe,试运行该语句，观察得出的结果，并利用help命令对你不熟悉的语句进行帮助信

息查

询，逐条给出上述程序段与结果的解释。

【求解】在MATLAB环境中感触如下语句，则可以看出，求解500£500随机矩阵的

逆，并

求出得出的逆矩阵与原矩阵的乘积，得出和单位矩阵的差，得出范数。•般来说，这样

得出

的逆矩阵精度可以达到10i12。

»tic,A=rand(500);B=inv(A);norm(A*B-eye(500)),toc3

ans=

1.2333e-012

Elapsedtimeis1.301000seconds.

2试用符号元素工具箱支持的方式表达多项式f(x)=x5+3x4+4x3+2x2+3x+

6,并令

x=sj1

s+1

,将f(x)替换成s的函数。

【求解】可以先定义出f函数，则由subs。函数将X替换成s的函数»symssx

f=x'5+3*x*4+4*x"3+2*x*2+3*x+6;

F=subs(f,x,(s-1)/(s+1))

F二

(s-l)5/(s+l)"5+3*(s-l)4/(s+1)4+4*(s-1)3/(s+1)3+

2*(s-l厂2/(s+l)-2+3*(s-l)/(s+l)+6

3用MATLAB语句输入矩阵A和B矩阵

①A二

2

664

1234

4321

23414

3241

3

775;

②

B

2

664

1+4j2+3j3+2j4+lj

4+lj3+2j2+3j1+4j

2+3j3+2j4+lj1+4j

3+2j2+3j4+lj1+4j

3

775

前面给出的是4£4矩阵，如果给出A(5;6)=5命令将得出什么结果？

【求解】用课程介绍的方法可以直接输入这两个矩阵

»A=[l234;4321;2341;3241]

A=

1234

6第2章MATLAB语言程序设计基础

43215

2341

3241

若给出A(5,6)=5命令，虽然这时的行和列数均大于B矩阵当前的维数，但仍然可以执

行该

语句，得出

»A(5,6)=5

A=

123400

432100

234100

324100

000005

复数矩阵也可以用直观的语句输入

»B=[l+4i2+3i3+2i4+1i;4+1i3+2i2+3il+4i;

2+3i3+2i4+lil+4i;3+2i2+3i4+lil+4i];

B二

1.0000+4.0000i2.0000+3.OOOOi3.0000+2.OOOOi4.0000+1.OOOOi4.0000+

l.OOOOi3.0000+2.OOOOi2.0000+3.OOOOi1.0000+4.OOOOi2.0000+3.OOOOi

3.0000+2.OOOOi4.0000+l.OOOOi1.0000+4.OOOOi3.0000+2.OOOOi2.0000+

3.OOOOi4.0000+l.OOOOi1.0000+4.OOOOi

4假设已知矩阵A,试给出相应的MATLAB命令，将其全部偶数行提取出来，赋给B矩

阵，6

用A二magic(8)命令生成A矩阵，用上述的命令检验一下结果是不是正确。

【求解】魔方矩阵可以采用magic()生成，子矩阵也可以提取出来»A=magic(8),

B=A(2:2:end,:)

A二

642361606757

955541213515016

1747462021434224

4026273736303133

3234352928383925

4123224445191848

4915145253111056

858595462631

B=

955541213515016

第2章MATLAB语言程序设计基础7

4026273736303133

4123224445191848

858595462631

5用MATLAB语言实现下面的分段函数y=f(x)二

8<

7

h;x>D

h=Dx;jxj6D

ih;x<jD

o

【求解】两种方法，其一，巧用比较表达式解决

»y=h*(x>D)+h/D*x.*(abs(x)<=D)-h*(x<-D);

另外一种方法，用循环语句和条件转移语句

»fori=l:length(x)

ifx(i)>D,y(i)=h;

elseifabs(x(i))<=D,y(i)=h/D*x(i);

else,y(i)=-h;end

end

其中，前者语句结构简单，但适用范围更广，允许使用矩阵型x,后者只能使用向量型

的

x,但不能处理矩阵问题。

6用数值方法可以求出S=

X63

i=0

2i=l+2+4+8+CCC+262+263,试不采用循环的形式求出和式的数值

解。由于数值方法采用double形式进行计算的，难以保证有效位数字，所以结果

不一定精确。试采用符号运算的方法求该和式的精确值。8【求解】用符号运算的方式

可以采用下面语句

»sum(sym(2)."[1:63])

ans=

18446744073709551614

由于结果有19位数值，所以用double型不能精确表示结果，该数据类型最多表示16

位有效

数字。其实用符号运算方式可以任意保留有效数字，例如可以求200项的和或1000项

的和可

以由下面语句立即得出。

>>sum(sym(2).[1:200])

ans=

3213876088517980551083924184682325205044405987565585670602750»

sum(sym(2).[1:1000])

ans=

214301721437253464189685009812000362112280962341106721488750077674070

210224987224498639675763139171625518934583510629365037429057138462808

8第2章MATLAB语言程序设计基础

719691551493971496078691355496484619708421492101247422837559083643060

9294996716388253479753511833108789215412582914239295537308433539208596

63305248773674411336138750

7编写一个矩阵相加函数matadd(),使其具体的调用格式为A二mat

add(A1,A2,A3,000),

要求该函数能接受任意多个矩阵进行解法运算。

【求解】可以编写下面的函数，用varargin变量来表示可变输入变量function

A=mat_acid(varargin)

A=0

fori=l:length(varargin),A=A+varargin{i};end

如果想得到合适的错误显示，则可以试用try,catch结构。function

A=matadd(varargin)

try

A=0;

fori=l:length(varargin),A=A+varargin{i};end

catch,error(lasterr);end

8自己编写一个MATLAB函数，使它能自动生成一个m£m的Hankel矩阵，并使其调

用格

式为v=[hl;h2;hm;hm+1;000;h2mj1];H=myhankel(v)o

【求解】解决这样的问题可以有多种方法：

①最直接的方法，Hi;j=hi+jjl,利用双重循环

functionH=myhankel(v)

(length(v)+l)/2;%严格说来还应该判定给定输入向量长度奇偶性lOfori=l:m,

forj=l:m

H(i,j)=v(i+j-l);

end,end

②考虑某一行(或列)，ai=[hi;hi+1;0M;hi+milL就可以用单重循环生成

Hankel矩阵了

functionH=myhankel(v)

m=(length(v)+l)/2;%严格说来还应该判定给定输入向量长度奇偶性for

H(i,:)=v(i:i+m-l);end

③利用现有的hankel()函数,则

functionH=myhanke1(v)

(length(v)+1)/2;%严格说来还应该判定给定输入向量长度奇偶性

H=hankel(v(l:m),v(m:end));

第2章MATLAB语言程序设计基础9

9已知Fibonacci数列由式ak=akj1+akj2;k=3;4;CCC可以生成，其中初值

为al=a2=1,

试编写出生成某项Fibonacci数值的MATLAB函数，要求

①函数格式为y=fib(k),给出k即能求出第k项ak并赋给y向量；②编写适当语

句，对输入输出变量进行检验，确保函数能正确调用；③利用递归调用的方式编写此函

数。

【求解】假设fib(n)可以求出Fibonacci数列的第n项，所以对>3则可以用

k=fib(ni

l)+fib(ni2)可以求出数列的n+1项，这可以使用递归调用的功能，11而递归调

用的出口为

U综上，可以编写出M-函数。

functiony=fib(n)

ifround(n)==n&n>=l

ifn>=3

y=fib(n-1)+fib(n-2);

else,y=1;end

else

error(，nmustbepositiveinteger.*)

encl

例如，n=10可以求出相应的项为

»fib(10)

ans=

55

现在需要比较一下递归实现的速度和循环实现的速度»tic,fib(20),toe

ans二

832040

elapsed_time=

62.0490

>>tic,a=[l1];fori=3:30,a(i)=a(i-l)+a(i-2);end,a(30),toeans=12

832040

elapsed_time二

0.0100

应该指出，递归的调用方式速度较慢，比循环语句慢很多，所以不是特别需要，解这样

问题

没有必要用递归调用的方式。

10第2章MATLAB语言程序设计基础

10由矩阵理论可知，如果一个矩阵M可以写成M=A+BCBT,并且其中A,B,C为

相应

阶数的矩阵，则M矩阵的逆矩阵可以由下面的算法求出

Mj1=

3

A+BCBT

，il

=Ai1iAiIB

3

Ci1+BTAjIB

il

BTAi1

试根据上面的算法用MATLAB语句编写一个函数对矩阵M进行求逆，并通过•个小例子

来

检验该程序，并和直接求逆方法进行精度上的比较。

13

【求解】编写这个函数

functionMinv=partinv(A,B,C)

Minv=inv(A)-inv(A)*B*inv(inv(C)+B'*inv(A)*B)*B'*inv(A);假设矩阵为

M=

2

664

51503616

50776032

36608748

16324868

3

775

且一知该矩阵可以分解成

A=

2

664

1000

0200

0030

0004

314

775

；B=

2

664

1234

2340

3400

4000

3

775

；C=

2

664

4000

0300

0020

0001

3

775

对这个例子。可以

»M=[51503616;50776032;36608748;16324868];iM=inv(M);%数值

逆，直接解法15

iM=

0.0553-0.03890.00170.0041

-0.03890.0555-0.0210-0.0021

0.0017-0.02100.0328-0.0137

0.0041-0.0021-0.01370.0244

»A=diag([l234]);B=hankel([l234]);C=diag([4321]);

iMl=part_inv(A,B,C)%分块矩阵的求解方法

iMl二

0.0553-0.03890.00170.0041

-0.03890.0555-0.0210-0.0021

0.0017-0.02100.0328-0.0137

0.0041-0.0021-0.01370.0244

乍看结果，似乎二者完全一致，实际上数值算法是有区别的。我们这里用解析方法得出

矩阵

的逆，然后用下面的语句比较两个结果的精度

第2章MATLAB语言程序设计基础11

>>Ml=sym(M);iM0=inv(Ml)

iMO=

[10713/193751,-7546/193751,332/193751,796/193751]

[-7546/193751,10759/193751,-4068/193751,-416/193751]

[332/193751,-4068/193751,19075/581253,-2652/193751]

[796/193751,-416/193751,-2652/193751,18919/775004]16

»norm(double(iMO)-iM)%直接求解的误差范数

ans=

2.7990e-017

>>norm(double(iMO)-iMl)%间接求解的误差范数

ans=

3.6583e-016

可见，用间接方法得出的逆矩阵误差更大，因为在这里新编写的函数中inv()函数使用

了多次，由此产生很大的传递误差。由此可以得出结论：如果某问题存在直接解，则尽

量别

使用间接方法，以加大传递误差。

11下面给出了一个迭代模型(

xk+1=1+yki1:4x2

k

yk+1=0:3xk

写出求解该模型的M-函数，如果取迭代初值为xO=0,yO=0,那么请进行30000次迭

代求出

一组x和y向量，然后在所有的xk和yk坐标处点亮一个点(注意不要连线)，最后绘

制出所

需的图形。提示这样绘制出的图形又称为Henon引力线图，它将迭代出来的随机点吸引

到一

起，最后得出貌似连贯的引力线图。

17

【求解】用循环形式解决此问题，可以得出如图2T所示的Henon引力线图。

»x=0;y=0;

fori=l:29999

x(i+l)=l+y(i)-l.4*x(i)-2;

y(i+l)=O.3*x(i);

end

plot(x,y,'.')

上述的算法由于动态定义x和y,所以每循环一步需要重新定维，这样做是很消耗时间

的，所以为加快速度，可以考虑预先定义这两个变量，如给出x二zeros(1,30000)。

12用MATLAB语言的基本语句显然可以立即绘制一个正三角形，试结合循环结构，编写

一个

小程序，在同一个坐标系下绘制出该正二角形绕其中心旋转后得出的一系列三角形，还

可以

调整旋转步距观察效果。

12第2章MATLAB语言程序设计基础

-1.5-1-0.500.511.5

-0.4

-0.3

-0.218

-0.1

0.1

0.2

0.3

0.4

图2THenon引力线图

【求解】假设正三角形逆时针旋转N度，则可以得出如图2-2a所示的示意图，三角形

的三个

顶点为（cos火sin四），（cos（N+120十）；sin（N+120十）），（cos（从+240十）；sin（以

+240十）），可以绘制出

其曲线，如图2-2b所示，试减小步距，如选择”=2;1;0:1,观察效果。

Nx

y

6

(a)示意图

-1-0.500.51

-1

-0.519

0.5

1

(b)曲线绘制效果

图2-2曲线绘制

»t=[0,120,240,0]*pi/180;%变换成弧度

xxx=[];yyy=[];

fori=0:5:360

tt=i*pi/180;

xxx=[xxx;cos(tt+t)];yyy=[yyy;sin(tt+t)];

end

第2章MATLAB语言程序设计基础13

plot(xxx',yyy','r'),axis('square*)

13选择合适的步距绘制出下面的图形sin

1

t

u

,其中t2(ii；Do

【求解】用普通的绘图形式，选择等间距，得出如图2-3a所示的曲线，其中x二0左

右显得

粗糙。

20

>>t=-l:0.03:1;y=sin(1./t);plot(t,y)

选择不等间距方法，可以得出如图2-3b所示的曲线。

»t=[-l:0.03:-0.25,-0.248:0.001:0.248,0.25:.03:1];y=sin(l./t);plot(t,y)

-1-0.500.51

-1

-0.5

0.5

1

(a)等间距曲线绘制

-1-0.500.51

-1

—0.5

0.5

1

(b)不等间距曲线绘制

图2-3不同自变量选取下的sin(l二t)曲线

14对合适的四范围选取分别绘制出卜列极坐标图形

①%=1:0013^2,②％=cos(71^=2),③％=sin(N)=%④为=1jcos3(7^)

【求解】绘制极坐标曲线的方法很简单，H]polar(^，^)即可以绘制出21极坐标图，

如图2-4所

示。注意绘制图形时的点运算：

»t=0:0.01:2*pi;subplot(221),polar(t,1.0013*t.2),%(a)

subplot(222),tl=0:0.01:4*pi;polar(tl,cos(7*tl/2))%(b)

subplot(223),polar(t,sin(t)./t)%(c)

subplot(224),polar(t,l-(cos(7*t)).*3)

15用图解的方式找到下面两个方程构成的联立方程的近似解。x2+y2=3xy2;x3j

x2=y2jy

【求解】这两个方程应该用隐式方程绘制函数ezplot。来绘制，交点即方程的解，如

图2-5a

所示。

14第2章MATLAB语言程序设计基础

20

40

30

210

60

240

90

270

120

30022

602402390

270

120

300

150

330

1800

1

2

30

210

60

240

90

270

120

300

150

330

1800

图2-4极坐标图

»ezplotCx*2+y*2-3*x*y2*);24holdon

ezplotCx-3-x2=y-2-y')

可用局部放大的方法求出更精确的值，如图2-5b所示。从图上可以精确读出两个交

点，(O:4O12;iO:8916),(1:5894;0:8185)。试将这两个点分别代入原始方程进行验

证。

—6—4—20246

-6

-4

-2

2

4

6

x

y

x3-x2=y2-y=0

(a)两个方程的曲线，交点为解

1.58941.58941.58941.58941.58941.58941.5894

0.8185

0.8185

0.818525

0.8185

0.8185

x

y

x3—x2=y2—y=0

(b)局部放大区域

图2-5二元联立方程的图解法

第2章MATLAB语言程序设计基础15

16请分别绘制出xy和sin(xy)的三维图和等高线。

【求解】(a)给出下面命令即可，得出的图形如图2-6a、b月示。»

[x,y]=meshgrid(-l1:1);

surf(x,y,x.*y),figure;contour(x,y,x.*y,30)

(b)给出下面命令即可，得出的图形如图2-6c、d所示。»[x,y]=meshgrid(-

pi1:pi);

surf(x,y,sin(x.*y)),figure;contour(x,y,sin(x.*y),30)

-1

1

-1

1

-126

-0.50

0.51

(a)xy三维图一1一0.500.51-1-0.50

0.51

(b)xy等高线一50

5

-50

5

-1-0.500.527

1

(c)sin(xy)三维图

-3-2-10123

-3

-2

-1

1

2

3

(d)sin(xy)等高线

图2-6三维图与等高线

17在图形绘制语句中，若函数值为不定式NaN,则相应的部分不绘制出来，试利用该规

律绘制

z=sinxy的表面图，并剪切下x2+y260:52的部分。

【求解】给出下面命令可以得出矩形区域的函数值，再找出x2+y260:52区域的坐

标，将其

函数值设辂成NaN,最终得出如图2-7所示的曲面。

16第2章MATLAB语言程序设计基础

>>[x,y]=meshgrid(-l1:1);z=sin(x.*y);

ii=find(x.2+y.^2<=0.52);z(ii)=NaN;surf(x,y,z)

-128

-0.5

0.5

1

-1

-0.5

0.5

1

-1

-0.5

0.5

1

图2-7得出的三维图

第3章微积分问题的计算机求解1试求出如下极限。①lim

x!l

(3x+9x)

1

x,②lim29

x!l

(x+2)x+2(x+3)x+3

(x+5)2x+5

【求解】极限问题由下面的语句可以直接求出。»symsx;f=(3求解〉］限/x);

limit(f,x,inf)ans=

9

»symsx;f=(x+2)-(x+2)*(x+3)-(x+3)/(x+5)-(2*x+5);limit(f,x,inf)

ans=

exp(-5)

2试求下面的双重极限。

①lim

x!j1

y!2

x2y+xy3

(x+y)3

,②lim

x!0

y!0

xy

p30

xy+1j1

,③lim

x!0

y!0

1icos

x2+y2

x2+y20

ex2+y2。

【求解】双重极限问题可以由下面语句直接求解。»symsxy

fa=(x2*y+x*y*3)/(x+y)*3;limit(1imit(fa,x,-1),y,2)ans=

-6

»fb=x*y/(sqrt(x*y+l)-l);limit(limit(fb,x,0),y,0)ans=

2

>>fc=(l-cos(x"2+y*2))*exp(x2+y*2)/(x"2+y*2);1imit(limit(fc,x,0),y,0)

ans=31

3求出下面函数的导数。

①y(x)=

q

xsinx

P

1iex,②y二

1i

P

cosax

x(1jcos

P

ax)

③atany

x

=ln(x2+y2),@y(x)=j

1

na

Inxn+a

xn;n>0

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

【求解】由求导函数diff()可以直接得出如下结果，其中③为隐函数，32故需要用

隐函数求

导公式得出导数。

»symsx;

f-sqrt(x*sin(x)*sqrt(1-exp(x)));simple(diff(f))

ans=

1/2/(x*sin(x)*(1-exp(x))(1/2))(1/2)*(sin(x)*(1-exp(x))(1/2)+x*cos(x)*(l-

exp(x))"(l/2)-l/2*x*sin(x)/(1-exp(x))*(l/2)*exp(x))»symsax

y=(1-sqrt(cos(a*x)))/(x*(l-cos(sqrt(a*x))))

simple(diff(y))

ans=

l/2/cos(a*x)(1/2)*sin(a*x)*a/x/(l-cos((a*x)(1/2)))-

(1-cos(a*x)"(l/2))/x2/(l-cos((a*x)"(l/2)))-l/2*(l-cos(a*x)"(1/2))/x/(1-

cos((a*x)"(1/2)))2*sin((a*x)"(1/2))/(a*x)(1/2)*a»f=atan(y/x)-log(x"2+y2);

fl=simple(-diff(f,x)/diff(f,y))

fl二

(y+2*x)/(x-2*y)

>>symsnpositive;symsa;f=-log((xn+a)/x"n)/(n*a);diff(f,x)ans=

一(n/x-(x-n+a)/(x-n)*n/x)/(x-n+a)*x八n/n/a

用LATEX表示上面的结果，则33

①1=2

从

sin(x)

P

1jex+xcos(x)p

1iexj1=2xsin(x)exp

1jex

11

1

q

xsin(x)

P

1iex

②1=2

sin(ax)ap

cos(ax)x(1jcos(p

ax))

1i34

P

cos(ax)x2(1jcos(pax))j1=231jp

cos(ax)'sin(pax)ax(1icos(pax))2pax③y+2xxj2y@j

35N

n

x

i

(xn+a)n

xnx

H

xn(xn+a)i1njlaj1

4试求出y(t)二

s

(xj1)(xj2)

(xj3)(xi4)

函数的4阶导数。

【求解】高阶导数可以由下面语句直接得出

»symsax

f=sqrt((x-1)*(x-2)/(x-3)/(x-4));simple(diff(f,x,4))

ans=

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

3*(16*x^ll-392*x"10+4312*x^9-28140*x^8+121344*x^7-

364560*x*6+783552*x'5-1214604*x*4+1342560*x*3-1015348*x*2+474596*x-

103741)/((x-l)*(x-2)/(x-3)/(x-4))*(7/2)/(x-3)*8/(x-4)~8336

16xllj392x10+4312x9；28140x8+121344x7i

364560x6+783552x5j1214604x4+1342560x311015348x2+474596xj103741』

(xi1)(xi2)

(xi3)(xj4)

W=2

(xi3)8(xj4)8

5在高等数学中，求解分子和分母均同时为0或1时，分式极限时可使用L'Mopital

法则，即

对分子分母分别求导数，再由比值得出，试用该法则lim

x!0

In(1+x)In(1jx)jln(ljx2)

x4

并和直接求出的极限结果相比较。

【求解】从给出的分母看，若想使之在X=0处的值不为0,则应该对其求4阶导数，

同样，

还应该对分子求4阶导数，将x=0代入结果，这样就可以使用L'Mopital法则求出

极限了。

>>symsx;n=log(l+x)*log(l-x)-log(l-x~2);d=x4;37

n4=diff(n,x,4);d4=diff(d,x,4);n4=subs(n4,x,0);L=n4/d4

L

1/12

现在直接求极限可以验证上述结果是正确的。»limit(n/d,x,0)

ans=

1/12

6已知参数方程

%

x=Incost

y=costjtsint

，试求出

dy

dx

和

d2y

dx2

t=p=3

o

【求解】参数方程的导数可以由下面语句直接求出。38»symst;x=log(cos(t));

y=cos(t)-t*sin(t);

diff(y,t)/diff(x,t)

ans=

-(-2*sin(t)-t*cos(t))/sin(t)*cos(t)

>>f=diff(y,t,2)/diff(x,t,2);subs(f,t,sym(pi)/3)

ans

3/8T/24*pi*3«l/2)

7假设u=cosi1

r

x

y

,试验证

@2u

@x@y

=@2u

@y@x

o

【求解】证明二者相等亦可以由二者之差为零来证明，故由下面的语句直接证明。

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

>>symsxy;u=acos(x/y);

diff(diff(u,x),y)-diff(diff(u,y),x)39

ans=

8设

(

xu+yv=0

yu+xv=1

,试求解

@2u

@x@y

【求解】用下面的语句可以直接得出如下结果。»symsxyuv

[u,v]=solve(，x*u+y*v=O','y*u+x*v=l','u,v');diff(diff(u,x),y)

ans=

2/(x-2-y八2厂2*x+8*y-2/(x-2-y-2厂3*x

9假设f(x;y)=

Zxy

ejt2dt,试求

x

y40

@2f

@x2

i2@2f

@x@y

+@2f

@y2

o

【求解】由下面的命令可以得出所需结果。»symsxyt

f=int(exp(-t2),t,0,x*y);

x/y*diff(f,x,2)-2*diff(diff(f,x),y)+diff(f,y,2)simple(ans)

ans二

-2*exp(-x-2*y-2)*(-x-2*y-2+1+x-3*y)

10假设已知函数矩阵f(x;y;z)

3x+eyz

x3+y2sinz

,试求出其Jacobi矩阵。

【求解】Jacobi矩阵可以由下面的语句直接得出。»symsxyz41

F=[3*x+exp(y)*z;x^3+y^2*sin(z)];jacobian(F,[x,y,z])

ans=

[3,exp(y)*z,exp(y)]

[3*x2,2*y*sin(z),y*2*cos(z)]

11试求解下面的不定积分问题。①I(x)=i

Z

3x2+a

x2(x2+a)2dx,②I(x)=Zp

x(x+1)

P

x+

P

1+x

dx

③I(x)二

Z

xeaxcosbxdx,@I(t)=Z

eaxsinbxsincxdx42

【求解】①该不定积分可以由下面的命令直接求出第3章微积分问题的计算机求解

21

>>symsxa

f=(3*x2+a)/(x2+(x^2+a)2);int(f,x)

ans二

12/(4+16*a)/(2+4*a+2*(1+4*a)(1/2))71/2)*

atan(2*x/(2+4*a+2*(l+4*a)Xl/2))'(1/2))+

48/(4+16*a)/(2+4*a+2*(1+4*a)71/2))71/2)*

atan(2*x/(2+4*a+2*(l+4*a)-(1/2)厂(1/2))*a+

12/(4+16*a)/(2+4*a+2*(1+4*a)X1/2)厂(1/2)*

atan(2*x/(2+4*a+2*(l+4*a)-(1/2))-(1/2))*(1+4*@厂(1/2)+

16/(4+16*a)/(2+4*a+2*(1+4*@厂(1/2)厂(1/2)*

atan(2*x/(2+4*a+2*(l+4*a)-(1/2))-(1/2))*(l+4*a)c(1/2)*a+12/(4+16*a)/(2+4*a-

2*(l+4*a)Xl/2))Xl/2)*

atan(2*x/(2+4*a-2*(l+4*a)-(1/2)厂(1/2))+

48/(4+16*a)/(2+4*a-2*(l+4*a)Xl/2)厂(1/2)*

atan(2*x/(2+4*a-2*(l+4*a)Xl/2))c(l/2))*a-

12/(4+16*a)/(2+4*a-2*(l+4*a「(1/2)厂(1/2)*

atan(2*x/(2+4*a-2*(l+4*a)-(1/2))-(l/2))*(l+4*a)-(1/2)-16/(4+16*a)/(2+4*a-

2*(l+4*a)Xl/2)),(l/2)*

atan(2*x/(2+4*a-2*(l+4*a)”1/2))Xl/2))*(l+4*a)-(l/2)*a②可以用下面的语句

求出问题的解43

>>symsx;f=sqrt(x*(x+1))/(sqrt(x)+sqrt(x+1));int(f,x);latex(ans)

并将其显示如下

2=15

P

x(x+l)x(3x+5)

P

x+1

i2=15

P

x(x+1)(x+1)(i2+3x)

P

x

③可以求出下面的结果

»symsabx

f=x*exp(a*x)*cos(b*x);int(f,x);latex(ans)其数学显示为

A

ax

a2+b2

a2ib244

(a2+b2)2

I

eaxcos(bx)j

A

bx

a2+b2+2ab

(a2+b2)2

!

eaxsin(bx)

④用下面的语句求解，得

>>symsxabc;f=exp(a*x)*sin(b*x)*sin(c*x);latex(int(f,x))

22第3章微积分问题的计算机求解亦即

1=2aeaxcos((bjc)x)

a2+(bic)2j1=2

(ib+c)eaxsin((bjc)x)

a2+(bjc)2

j1=2aeaxcos((b+c)x)

a2+(b+c)2+1=2

(ibjc)eaxsin((b+c)x)45

a2+(b+c)2

12试求出下面的定积分或无穷积分。

①I=

Z1

cosX

P

X

dx,②I二

Z1

1+x2

1+x4dx

【求解】①可以直接求解

>>symsx;int(cos(x)/sqrt(x),x,0,inf)

ans=

1/2*2、(1/2)*pi-(l/2)

②可以得出

>>symsx;int((l+x*2)/(l+x*4),x,0,1)

ans=

1/4*271/2)*pi

13假设f(x)=ej5xsin(3x+p=3),试求出积分函数R(t)=46Zt

f(x)f(t+x)dxo

【求解】定义了x的函数，则可以由subs()函数定义出t+x的函数，这样由下面的

语句可

以直接得出R函数。

»symsxt;f=exp(-5*x)*sin(3*x+sym(pi)/3);

R=int(f*subs(f,x,t+x),x,0,t);simple(R)

ans二

1/1360*(15*exp(t)TO*3c(1/2)Mos(3*t)-25*cos(9*t)+

25*exp(t)"10*3"(1/2)*sin(3*t)-68*cos(3*t)-15*3"(1/2)*cos(9*t)-

25*3X1/2)*sin(9*t)-15*exp(t)10*sin(3*t)+15*sin(9*t)+

93*exp(t)10*cos(3*t))/exp(t)15

该结果可以写成

1

1360

15(et)10p

3cos(3t)j68cos(3t)j15(et)10sin(3t)j25

P

3sin(9t)+25(et)10p

3sin(3t)

+15sin(9t)j25cos(9t)i1547

P

3cos(9t)+93(et)10cos(3t)

(et)15

14对a的不同取值试求出I=

Z1

cosax

1+x2dxo

【求解】由下面的循环结构可以得出不同a值下的无穷积分值，并可以绘制出它们之间

关系

的曲线，如图3T所示。

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

>>symsxa;f=cos(a*x)/(l+x2);aa=[0:0.l:pi];y=[];

forn=aa

b=int(subs(f,a,n),x,0,inf);y=[y,double(b)];

end

plot(aa,y)

00.511.522.533.5

0.2

0.4

0.648

0.8

1

1.2

1.4

1.6

图3T不同a值下的积分值曲线

15试对下面函数进行Fourier幕级数展开。

①f(x)=(pijxj)sinx;ip6x<p②f(x)=ejxj;ip6x<p③f(x)=

(

2x=l;0<x<1=2

2(1jx)=l;1=2<x<1

,且1二p。

【求解】①可以立即由下面的语句求出。

>>symsx;f=(sym(pi)-abs(x))*sin(x);

[A,B,F]=fseries(f,x,10,-pi,pi);F

F二

l/2*pi*sin(x)+16/9/pi*sin(2*x)+32/225/pi*sin(4*x)+

48/1225/pi*sin(6*x)+64/3969/pi*sin(8*x)+80/9801/pi*sin(10*x)该结果在LATEX

下可以显示为

1

2psinx+49

0.8

1

1.2

1.4

1.6

图3T不同a值下的积分值曲线

15试对下面函数进行Fourier基级数展开。

①f(x)=(pjjxj)sinx;ip6x<p②f(x)=ejxj;ip6x<p③f(x)=

(

2x=l;0<x<1=2

2(1jx)=l;1=2<x<1

.且1=p。

【求解】①可以立即由下面的语句求出。

»symsx;f=(sym(pi)-abs(x))*sin(x);

[A,B,F]=fseries(f,x,10,-pi,pi);F

F=

l/2*pi*sin(x)+16/9/pi*sin(2*x)+32/225/pi*sin(4*x)+

48/1225/pi*sin(6*x)+64/3969/pi*sin(8*x)+80/9801/pi*sin(10*x)该结果在LATEX

下可以显示为

1

2psinx+49

sinlOx

P

②可以由下面语句求解，并得出数学公式为»symsx;f=exp(abs(x));

[A,B,F]=fseries(f,x,10,-pi,pi);F

得出的解析解为

1=2

2epj2

P

+

(iepj1)cos(x)

P

+

(2=5epi2=5)cos(2x)

P

+

(j1=5epj1=5)cos(3x)

P

24第3章微积分问题的计算机求解+

(2=17epi2=17)cos(4x)

p51

(j1=13epi1=13)cos(5x)p

+

i2

37epj2

37

0

cos(6x)

P

+

(i1=25epi1=25)cos(7x)p

+

i2

65epj2

65

0

cos(8x)

P

+

(j1=41epi1=41)cos(9x)52p

+

i2

101epj2

101

cos(10x)

P

进一步观察结果可见，该式子可以手工化简，例如提取系数(epi

(iDn)=p。或对各项系数

逐项求值(保留10位有效数字)

»vpa(F,10)

ans=

7.047601355-7.684221126*cos(x)+2.819040541*cos⑵*x)T.536844225*cos(

3.*x)

+.8291295709*cos(4.*x)5910939328*cos(5.*x)+.3809514246*cos(6.*x)-

.3073688450*cos(7.*x)+.2168492724*cos(8.*x)-.1874200274*cos(9.*x)

+.1395564625*cos(10.*x)

③似乎求解起来更困难，巧妙利用符号运算工具箱中的heavisideO函数，则可以将

原函数

表示成

f(x)=20heaviside53

3

xip

2

2

P

x

JXip=2j

xip=2

这样就可以用下面的语句求出函数的Fourier级数。

»symsx;pil=sym(pi);

f=2*heaviside(x-pi1/2)-2/pil*x*abs(x-pil/2)/(x-pil/2);

[a,b,F]=fseries(f,x,10,-pi,pi);F

F=

-1/4+4/pi2*cos(x)+(4/pi+2)/pi*sin(x)-2/pi"2*cos(2*x)-l/pi*sin(2*x)+

4/9/pi-2*cos(3*x)+(-4/9/pi+2/3)/pi*sin(3*x)T/2/pi*sin(4*x)+

4/25/pi/2*cos(5*x)+(4/25/pi+2/5)/pi*sin(5*x)-2/9/p12*cos(6*x)-

l/3/pi*sin(6*x)+4/49/pV2Mos(7*x)+(-4/49/pi+2/7)/pi*sin(7*x)-

l/4/pi*sin(8*x)+4/81/pi八2Mos(9*x)+(4/81/pi+2/9)/pi*sin(9*x)-

2/25/pi^2*cos(10*x)-l/5/pi*sin(10*x)

i1=4+454

cos(x)p2+j

4pj1+20

sin(x)p

i2

cos(2x)p2j

sin(2x)p

+4=9cos(3x)p2+j

i4=9pi1+2=30

sin(3x)p55

i1=2sin(4x)p+425

cos(5x)p2+j4

25pi1+2=50

sin(5x)pj2=9cos(6x)p2j1=3sin(6x)p+44956cos(7x)p2+jj4

49pi1+2=70

sin(7x)pi1=4sin(8x)p+481

cos(9x)p2+j4

81pj1+2=90

sin(9x)pj572

25

cos(lOx)

P2

i1=5

sin(lOx)

P

16试求出下面函数的Taylor基级数展开。第3章微积分问题的计算机求解25①

Zx

sint

t

dt②In

N

1+x

1ix

③In

3

x+58

P

1+x2

④(1+4:2x2)0:2

⑤ej5xsin(3x+p=3)分别关于x=0>x=a的累级数展开⑥对f(x;y)=

1jcos

i

x2+y2

0

i

x2+y2c

ex2+y2关于x=1;y=0进行二维Taylor暴级数展开。

【求解】由下面的语句可以分别求出各个函数的嘉级数展开，由latex(ans)函数可以

得出下

面的数学表示形式。

>>symstx;f=int(sin(t)/t,t,0,x);

taylor(f,x,15)

>>symsx;f=log((l+x)/(l-x))»taylor(f,x,15)

>>symsx;f=log(x+sqrt(1+x2));taylor(f,x,15)

>>symsx;f=(l+4.2*x2)0.2;taylor(f,x,13)

①xi1=18x3+59

1

600x5i

1

35280x7+

1

3265920x9i

1

439084800x11+

1

80951270400x13

②2x+2=3x3+2=5x5+2=7x7+2=9x9+2=11x11+2=13x13③xj1=6x3+

3

40x5i

5

112x7+

35

1152x9i

63

2816x11+

231

13312x1360

④1+

21

25x2i

882

625x4+

55566

15625x6j

4084101

390625x8+

1629556299

48828125xlOi

136882729116

1220703125xl2

⑤该函数的前4项展开为

»symsxa;f=exp(-5*x)*sin(3*x+sym(pi)/3);taylor(f,x,4,a)

ej5asin

3

3a+p

3

+61

3

3ei5acos3

3a+p3'j5ej5asin3

3a+p3

(xia)+3

8ej5asin3

3a+p3'j15ej5acos33a+p62

3

(xia)2

+

人人文库> 全部分类> 教育资料 > 辅导培训

温馨提示

  • 1. 本站所有资源如无特殊说明，都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等，如果需要附件，请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸，网页内容里面会有图纸预览，若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对用户上传分享的文档内容本身不做任何修改或编辑，并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容，请与我们联系，我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

最新文档

评论

0/150

提交评论