MATLAB在高等数学中的应用

MATLAB 在高等数学中的应用
１ 引言
在科学技术不断发展的今天，计算机得到迅速发展．计算机的出现归功于数学家的奠基性工作，计算机的发展又为数学的发展提供了威力无比的武器和工具，从而彻底改变了长期以来数学仅仅靠一支笔，一张纸的传统，使数学的应用在广度及深度两方面都达到了前所未有的程度，深刻地影响了数学的发展进程和思维模式，同时也使数学技术成为现今高科技的一个重要组成部分和突出标志． 中国科学院院士王梓坤在《今日数学及其应用》一文中指出“精确定量思维是对21世纪科技人员共同的素质要求。

所谓定量思维就是指人们从实际问题中提炼数学问题，抽象化为数学模型，用数学计算求出此模型的解或近似解，然后回到现实中进行检验，必要时修改模型使之更切合实际，最后编制解决问题的软件包，以便得到更广泛的方便的应用”。

在当前众多数学应用软件中，MATLAB 是一个应用广泛、功能强大的软件．在７０年代后期，Cleve Morler 博士编写了MATLAB ．１９８４年， Cleve Morler 和John Little 成立Math Works 公司，正式把MA TLAB 推向市场，并对MA TLAB 进行不断的开发．MATLAB 已经发展成为适合多学科的功能强大的大型软件．在欧美等高校，MA TLAB 已经成为线性代数，自动控制理论，数理统计，数字信号处理，动态仿真等高级课程的基本教学工具，同时被研究单位和工业部门广泛应用，使科学研究和解决各种具体问题的效率大大提高．MA TLAB 提供了专业水平的数值计算，符号计算和图形可视化等功能，它几乎可以解决实际应用中出现的绝大多数的数值计算问题，如数据分析，曲线拟合，数值分析等．MA TLAB 软件不仅能够进行简单的数值计算，还能进行求导，积分，解方程，求特征值和特征向量等符号计算，并且MA TLAB 的图形功能强大，既包括对二维和三维数据可视化，图像处理，动画制作等高层次的绘图命令，也包括可以完全修改图形局部及编制完整图形界面的低层次的绘图命令．
MATLAB 作为数学软件用于解决高等数学中一些计算问题和绘图问题，给学生一种全新的感觉，激发起学习的兴趣，加深对所学知识的理解，使学生对数学发展现状及应用有切实的体会． 如在高等代数中，矩阵的幂方和除法是两个计算量比较大而且容易出错的运算，尤其是幂方，而这些在MATLAB 中都会很快又准确的得出结果．
例１ 已知A=⎪⎪⎪⎭
⎫ ⎝⎛01410141061062，求A 3．
应用MATLAB 命令：
A=[2 6 10;6 10 14;10 14 0]^3
ans=
2592 4416 4368
4416 7392 6768
4368 6768 3840
上面的例题已经体现出MATLAB 的优势．同样在令很多人头疼的绘图上MA TLAB 也能很快给出符合要求而精确的图形．
例２ 作出螺旋线t z t y t x ===,cos ,sin 的图形．
应用MATLAB 命令：
>> t=0:pi/50:10*pi; >>plot3(sin(t),cos(t),t)
通过上面的例题使学生在心理上对高等数学的学习不再感到枯燥，将书本上学习的高等数学知识变成计算机语言很快就能得出结果，即使是需要费时的计算和绘图，这些都有助于学生自觉的将MATLAB 应用在高等数学的学习中．以下介绍MA TLAB 在高等数学中的应用．
２ 微积分方面
数值分析主要讨论有关函数的极值问题，零点问题以及微积分问题等，MATLAB 提供了大量的解决上述数值分析任务的函数．
２.１ 极值函数
MATLAB 的主要极值函数是fminbnd ，该函数的常用调用格式为：
X=fimnbnd(FUN,Ｘ１,Ｘ２)
其中输入参数FUN 为欲求极值的函数，Ｘ１，Ｘ２为求极值的范围．
例３求函数x
cot在［１，２］上的极值．
应用MATLAB命令：
>>y=fminbnd('cot',1,2)
>>y=
1.9999
如何求函数的极大值？思考一下．
２.２求根函数
函数fzero为求根函数，其调用格式为：
X=fzero(FUN,X0)
其中输入参数FUN为欲求根的函数，当X0为一个数值时，将得到X０最近的方程的根，当X0为两个数值时，将求出两个元素之间的方程的根．
例４求函数x
cot在１附近的零点，在区间［１，２］中的零点．
应用MATLAB命令：
＞＞x=fzero('cot',1)
＞＞x=
1.5708
＞＞x=fzero('cot',[1,2])
＞＞x=
1.5708
２.３数值积分函数
MATLAB中的数值积分函数是quad,其调用格式为：
Q=quad(FUN,A,B)将求出函数FUN在A和B之间的积分结果．
例５求
2
1
sin xdx
．
应用MATLAB命令：
>> q=quad('sin',1,2)
q =
0.9564
２.４微分函数
Diff(s)对符号表达式s对自变量的积分．
Diff(s,'v')对以v为自变量的符号表达式s求微分．
Diff(s,'v',n) 其中n 为正整数，对函数表达式s 求n 阶微分．
例６ 求sint 2的一阶微分，t 7的七阶微分．
应用MATLAB 命令：
>> t=sym('t');
>> u=sym('u');
>>diff(sin(t^2))
ans =
2*cos(t^2)*t
>> diff(t^7,7)
ans =
5040
在高等数学中函数列的一致收敛性概念比较抽象难理解，而借助MATLAB 强大的绘图功能就能形象直观的理解一致收敛性并对函数列的一致收敛性进行判别．
例７ 已知如下三个函数序列：
];5,0[},{)}({2∈=-x e nx x f nx n
];1,0[},{)}({2
22∈=-x xe n x g x n n
).,0[},)
1(1{)}({2222
2+∞∈+-=x x n x n x h n 通过作出它们的图象，根据一致收敛的几何意义来观察它们的一致收敛性 .
(1) 作出)}({x f n 的一族曲线如下：
应用MATLAB 命令：
hold on; x=0:0.02:5; for n=1:8
y=n*x.^2.*exp(-n.*x); plot(x,y); end;
图1
由图1可见, )}({x f n 中每条曲线的峰值随着∞→n 而趋于0 , 故)}({x f n 能在[ 0, 5 ]上一致收敛于0．
(2) 作出)}({x g n 的一族曲线如下：
应用MATLAB 命令：
clf;
hold on;
x=0:0.01:1;
for n=5:-1:1
y=n^2*x.*exp(-n^2*x.^2);
plot(x,y);
end
图2 由图2可见，
)}({x g n 中每条曲线的峰值随着∞→n 而趋于∞+ 故在0=x 的任意小右邻域内, )}({x g n 不可能一致收敛；但是对任意小的正数a , )}({x g n 在[a , 1 ]上仍有可能一致收敛于0．
(3) 作出)}({x h n 的一族曲线如下：
应用MATLAB 命令：
>>clf;
>>hold on;
>>x=0:0.01:3;
>>for n=1:1:5
y=(1-n^2*x.^2)./(1+n^2*x.^2).^2;
plot(x,y);
end ;
图3
由图3可见，极限函数在0=x 处不连续，故在0=x 的任意小右邻域内, )}({x h n 不可能一致收敛；但在任何),0[),[+∞⊂+∞a 上，)}({x h n 仍有可能一致收敛于0．
在高等数学中幂级数求和函数手算比较困难，应用MATLAB 求和函数就会很快求出和函数．
例８ 求幂级数∑∞
=-2)1(k k
k k x 的和函数S2．
应用MATLAB 命令：
syms k x
S2=symsum('x^k/k/(k-1)',k,2,inf)
S2 = 1/2*x^2*(2/x*(-log(1-x)/x-1)
*(x-1)-1/(x-1)*(-2*x+2))
S2=simplify(S2) (将S2简化)
S2 = -log(1-x)*x+log(1-x)+x
在上面的例题中应用MATLAB 对高等数学中出现的一些复杂计算和函数项序列一致收敛进行判别，体现了MATLAB 在高等数学中的应用．
３ 矩阵
MATLAB 的数学能力大部分是从它的矩阵函数中派生出来的，在MA TLAB 中提供了丰富的矩阵运算函数，包括所有的标准的线性代数中使用到的矩阵函数．
３.１ 计算行列式、秩
Det 计算矩阵行列式的值
Rank 函数计算矩阵的秩． n=1
n=5
例９ A=⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛1097591086781075675,求｜Ａ｜，rank(A). 应用MATLAB 命令：
>> A=[5 7 6 5 ;7 10 8 7;6 8 10 9;5 7 9 10]；
>> det(A)
ans =
1
>> rank(A)
ans =
4
３.２ 特征值的计算
eig （Ａ） 计算矩阵的特征值
如果Ａ是n n *矩阵，若λ满足ＡＸ＝λＸ，则称λ为Ａ的特征值．
例10 Ａ＝ ⎪⎪⎪⎪⎪⎭
⎫ ⎝
⎛86109638715435442 ,求A 的特征值. 应用MATLAB 命令： >>A=[2 4 4 5;3 4 5 1;7 8 3 6;9 10 6 8]；
＞＞eig=(A)
ans =
20.4983
-3.4064
-1.3133
1.2214
如果还要求特征向量，则可以用eig （Ａ）函数的第一个返回值得到．
应用MATLAB 命令：
>> [X,D]=eig(A)
X =
-0.3659 -0.6486 0.7050 0.4314
-0.2674 -0.2135 -0.6239 -0.7569
-0.5175 0.6321 0.2781 0.0644
-0.7258 0.3664 -0.1905 0.4868
D =
20.4983 0 0 0
0 -3.4064 0 0
0 0 -1.3133 0
0 0 0 1.2214
在MATLAB 中特征值和特征向量的求解过程体现了矩阵作为MA TLAB 的基本计算单位在进行矩阵运算时是很简便的．
４ 方程和方程组
“方程是很多工程和科学工作的发动机”．若干世纪以来，工程师和科学家们花了大量的时间用于求解方程（组）．数学家研究各种各样的方程求解方法．MATLAB 为求解方程（组）提供了便利条件．
４.１ 方程
MA TLAB 的命令输入格式：
solve('eqn1','eqn2'…,'eqnN','var1','var2',…'varN')
其中eqni 表示第i 个方程，vari i 表示第i 个变量，i=1,2,…,N ．
在求解多项式方程时还可以用MATLAB 函数roots(p)，其中p 是多项式的系数按降幂排列所形成的n+1维列向量，它能够给出全部根（包括重根）．避免在高等数学日常计算中出现的丢根情况．
例11 求解多项式方程018
9=++x x .
应用MATLAB 命令：
>> P=[1,1,0,0,0,0,0,0,0,1];
>> roots(P)
ans =
-1.2131
-0.9017 + 0.5753i
-0.9017 - 0.5753i
-0.2694 + 0.9406i
-0.2694 - 0.9406i
0.4168 + 0.8419i
0.4168 - 0.8419i
0.8608 + 0.3344i
0.8608 - 0.3344i
４.2 方程组
４.2.1 线性方程组
除了使用MA TLAB 函数solve 以外，还可以用其他MATLAB 命令．如将线性方程组写成矩阵形式A*X=b ，其中Ａ为m*n 阶方阵，X 和b 均为n 阶列向量．
例12 解线性方程组⎪⎪⎪⎪⎪⎭⎫ ⎝
⎛852*********X=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛514351804366 .
应用MATLAB 命令： >> A=[1 2 3;4 5 6;7 8 0;2 5 8];
>>b=[366 804 351 514]';
>> x=A\b
x =
247.9818
-173.1091
114.9273
４.2.2 非线性方程组
非线性方程组可以用solve( )函数进行求解，一般给出的是数值解，也可以用fsolve( )进行求解，输入格式为
fsolve('FUN',X0)
其中FUN 表示M 文件函数，X0表示变量的初始点．下面用MATLAB 求解一个比较复杂的非线性方程组．
例13 解非线性方程组⎪⎩
⎪⎨⎧=-++=+-+=-++050
12307ln sin 32z y x z x z y x y .
应用MATLAB 命令：
function eq=nxxf(x)
global number;
number=number+1;
eq(1)=sin(x(1))+x(2)^2+log(x(3))-7;
eq(2)=3*x(1)+2^x(2)-x(3)^3+1;
eq(3)=x(1)+x(2)+x(3)-5;
global number;
number=0;
y=fsolve('nxxf',[1,1,1])
y =
0.5991 2.3959 2.0050
５ 微分方程
微分方程是高等数学中一个比较重要的分支,它的应用也比较广泛,使用MATLAB 软件就可以方便地求出它的各种形式的解．
５.１ 解析解
求微分方程的解析解的MATLAB 命令：
dsolve('eqn1','eqn2',…,'x')
其中'eqni'表示第i 个方程，'x'表示微分方程组中的自变量．
除了用MA TLAB 求解一阶微分方程外，也可以解决较复杂的二阶微分方程．
例14 求解二阶微分方程
2x y ''x +y '2/1,/2)2/(),2/(,0)(22==-'==-+n y y y n x πππ．
应用MA TLAB 命令：
>> dsolve('D2y+(1/x)*Dy+(1-(1/2)^2/x^2)*y=0','y(pi/2)=2','Dy(pi/2)=(-2/pi)','x')
ans =
2^(1/2)*pi^(1/2)/x^(1/2)*sin(x)
５.２ 数值解
在没有解析解的情况下就要求数值解．
求解微分方程（组）数值解，用２阶龙格－库塔公式和４阶龙格－库塔公式，其程序为：
[t,y]=ode23('F',ts,y0,options) [t,y]=ode45('F',ts,y0,options)
其中F 为由微分方程（组）写成的M 文件名．
例15 求方程0)(2
22=-+'+''y n x y x y x 的数值解.
应用MATLAB 命令： function f=jie3(x,y)
f=[y(2);-y(2)/x+((1/2)^2/x^2-1)*y(1)]; >> [x,y]=ode23('jie3',[pi/2,pi],[2,-2/pi]) x = 1.5708 1.6074 1.7645 1.9215 2.0786 2.2357 2.3928 2.5499 2.7069 2.8640 3.0211 3.1416 y =
2.0000 -0.6366 1.9758 -0.6869 1.8518 -0.8879 1.6982 -1.0631 1.5192 -1.2108 1.3193 -1.3293 1.1032 -1.4174 0.8756 -1.4744 0.6416 -1.5002 0.4060 -1.4951 0.1735 -1.4602
0.0002 -1.4140
>> y1=y(:,1);
>> y2=y(:,2);
>> plot(x,y1,x,y2,'r'),gtext('y1'),gtext('y2')
６插值
MATLAB中的插值函数为interpl( ),其调用格式为
yi=interpl(x,y,xi,'method')
其中x,y为插值点，yi为在被插值点xi处的插值结果；x,y为向量．'method'表示采用的插值方法．MA TLAB提供的插值方法有几种：'nearest'最邻近插值;'linear'线性插值;'spline'三次样条插值;'cubic'立方插值．用MATLAB插值做一个数学实验．
O）例16 在一天24小时内，从零点开始每间隔２小时测得的环境温度数据为（C
12，9，9，1，0，18，24，28，27，25，20，18，15，13，
推测中午１点（即13点）时的温度．
应用MATLAB命令：
>> x=0:2:24;
>> y=[12 9 9 10 18 24 28 27 25 20 18 15 13];
>> x1=13;
>> y1=interp1(x,y,x1,'spline')
y1 =
27.8725
O
所以，中午１点时的温度大约为27.8725C
若要得到一天24小时的温度曲线，只需继续键入：
应用MATLAB命令：
>> xi=0:1/3600:24;
>> yi=interp1(x,y,xi,'spline');
>> plot(x,y,'o',xi,yi)
７曲线拟合
MATLAB中的拟合函数为polyfit( ),其调用格式为:
Polyfit(X,Y,n)
执行该函数将产生一个n阶多项式P，并且使得P(X)=Y.
例17下表展示了来自１９８５年以来第x年从中国进口到美国的（以百万美元计数）值y （a）求一能近似该数据的线性函数y＝b＋mx．
（b）对每个y求lny，并利用y和lny的值求一个能近似该数据的线性函数lny＝b＋mx．
e的指数函数．（c）利用（b）中的方程求一个能拟合该数据的形为y＝a kx
年1985 1986 1987 1988 1989 1990 1991
第x年0 1 2 3 4 5 6
进口值y 3862 4771 6293 8511 11990 15237 18976
应用MATLAB命令:
(a)线性拟合
x=[0:6];
>> y=[3862 4771 6293 8511 11990 15237 18976];
>> plot(x,y,'o')
>> hold on
>> f1=polyfit(x,y,1)
f1 =
1.0e+003 *
2.5704 2.2374 >> y1=polyval(f1,x); >> plot(x,y1,'r')
线性拟合函数拟合效果
(b) 对数线性拟合 f2=polyfit(x,log(y),1) f2 =
0.2765 8.2275 >> y2=polyval(f2,x); >> hold on >> plot(x,log(y),'o') >> plot(x,y2,'k')
对数线性函数拟合效果
(c) 由(b)可知，x
e
y 2756.05.3742
plot(x,y,'o')
>> hold on
>> plot(x,3742.5*exp(0.2756*x))
指数拟合效果
由拟合效果知指数函数拟合较好．
８数学实验
数学实验是对实际生活中的数学问题的解决，提出问题并进行探究、思考、分析等思维活动，最终理解并解决问题，提出想法借助于数学软件实现．培养学生的独立能力和创新精神．波音公司飞机最佳定价策略[2](P46-48)
全球最大的飞机制造商—波音公司自1955年推出波音707开始．成功地开发了一系列的喷气式客机．问题：讨论该公司对一种新型客机最优定价策略的数学模型．
(1) 问题分析
定价策略涉及到诸多因素，这里考虑以下主要因素：
价格、竞争对手的行为、出售客机的数量、波音公司的客机制造量、制造成本、波音公司的市场占有率等等因素。

(2) 假设及模型
价格记为p，根据实际情况，对于民航飞机制造商，能够与波音公司抗衡的竞争对手只有一个，因此他们可以在价格上达成一致，具体假设如下：
1）型号：为了研究方便，假设只有一种型号飞机；
销售量：其销售量只受飞机价格p的影响．预测以此价格出售，该型号飞机全球销售量为Ｎ．Ｎ应该受到诸多因素的影响，假设其中价格是最主要的因素．根据市场历史的销售规律和需求曲线，假设该公司销售部门预测得到N=N(p)=-78p2+655p+125.
(3) 市场占有率：既然在价格上达成一致，即价格的变化是同步的，因此，不同定价不会影响波音公司的市场占有率，因此市场占有率是常数，记为h ．
4) 制造数量：假设制造量等于销售量，记为x ．既然可以预测该型号飞机全球销售量，结合波音公司的市场占有率，可以得到 x=h*N(p).
5) 制造成本：根据波音公司产品分析部门的估计，制造成本为 C(x)=50+1.5x+8x 4
3. 6) 利润：假设利润等于销售收入去掉成本，并且公司的最优策略原则为利润R （p)最大．利润函数为 R(p)=px-C(x).
由以上简化的分析及假设得到波音公司最佳定价策略的数学模型如下 MaxR(p)=px-C(x),
其中N(p)=-78p 2
+655p+125,x=h*N(p),C(x)=50+1.5x+8x 4
3
,p,x,N ≥0. （３）模型求解
我们采用图形的方法求解．最优定价策略下价格p 大致在6到7之间．p ≈6.2859,R=1780.8336. MATLAB 程序如下:
>> h=0.5;a=0;b=8;n=80;d=(b-a)/n; >>for i=1:n+1 pr(i)=a+(i-1)*d;p=pr(i); >>n=-78*p^2+655*p+125;x=h*n;r=p*x; >>c=50+1.5*x+8*x^(3/4);l(i)=r-c; >>end
>>plot(pr,l);grid on
>>xlabel('价格p');title('利润曲线R(p)')
通过本文体现了Matlab在高等数学中的应用是很广泛的，对高等数学的学习起到了很好的辅助作用．
参考文献：
[1] 云舟工作室.MATLAB 6 数学建模基础教程[M].北京:人民邮电出版社,2001
[2] 付鹏.数学实验.北京:科学出版社,2000
[3] Eric Connally Hughes-Hallett et al.,Functions modeling change:Apreparation for calculus,John Wiley&Sons,InC.,2000;Preliminary edition,1998.
[4] 陈桂明.应用MATLAB建模与仿真[M].北京:科学出版社,2001
[5] 陈桂明.MATLAB数理统计(6.X)[M].北京:科学出版社,2002
[6] 苏金明.MATLAB 6.1实用指南[M].北京:电子工业出版社,2002
[7]G.Fulford,Modelling with differential and difference equations,Cambridge: Cambridge University press,1997.。