《数学实验与Matlab》程序.docx
《《数学实验与Matlab》程序.docx》由会员分享，可在线阅读，更多相关《《数学实验与Matlab》程序.docx（91页珍藏版）》请在冰豆网上搜索。
《数学实验与Matlab》程序
目录
《数学实验与Matlab》程序1
实验１．矩阵运算与Matlab命令2
1.1知识要点与背景：
日常矩阵及其运算2
1.2实验与观察：
矩阵和Matlab语言2
1.2.1向量的生成和运算2
1.2.2.矩阵创建和运算3
1.2.3Matlab工作环境和编程5
1.3应用、思考与练习5
1.3.1关系矩阵6
1.3.2投入产出6
1.3.3循环比赛的名次6
1.3.4参考程序6
实验２．函数的可视化与Matlab作7
2.1实验与观察：
函数的可视化7
2.1.1Matlab二维绘图命令7
2.1.2多元函数的可视化与空间解析几何（三维图形）10
2.2应用、思考和练习12
2.2.1线性p周期函数12
2.2.2平面截割法和曲面交线的绘制12
2.2.3微分方程的斜率场13
2.2.4颜色控制和渲染及特殊绘图指令14
实验３．函数式－直接确定型模型15
3.1知识要点和背景：
函数—直接确定性模型15
3.2实验与观察：
插值与拟合15
3.2.1插值方法与多项式拟合的概念16
3.2.2用Matlab作插值和拟合16
3.2.3用鼠标选节点观察插值、拟合的效果16
3.2.4观察程序说明16
3.3应用、思考和练习17
3.3.1若干函数的插值和拟合练习18
3.3.2几个应用问题18
3.3.4多元函数的插值18
3.3.5Fourier级数和周期函数的经验公式18
3.3.6由实验到理论:
从开普勒定律和牛顿万有引力定律19
实验４.微分、积分和微分方程19
4.1.知识要点和背景：
微积分学基本定理19
4.2实验与观察（Ⅰ）：
数值微积分19
4.2.1实验：
积分定义、微分方程和微积分基本定理的联系19
4.2.2求解数值积分的Matlab积分命令20
4.2.3观察程序及其说明20
4.3实验与观察（Ⅱ）：
Matlab符号微积分简介22
4.3.1创建符号变量22
4.3.2求符号极限limit指令22
4.3.3求导指令diff22
4.3.4求符号积分int23
4.3.5化简、提取和代入23
4.4应用、思考和练习24
4.4.1追击问题24
4.4.2应用问题24
实验５.用导数作定性分析25
5.1知识要点：
函数作图—用导数定性描述函数25
5.2实验与观察：
微分方程的定性解图示26
5.2.1人口增长的预测26
5.3应用、思考和练习28
5.3.1.函数作图29
5.3.2.平衡点的分类29
5.3.3定性分析的应用29
实验６．最优化实验30
6.1知识要点与背景30
6.1.1由简入繁：
最佳水槽断面问题的推广30
6.1.2微分法求最大和最小30
6.2实验与观察（Ⅰ）：
模拟盲人下山的迭代寻优法31
6.3.实验与观察（Ⅱ）：
Matlab优化工具箱简介32
6.3.1多元函数无约束优化指令fminunc和fminsearch32
6.3.2其它的优化算法指令37
6.4应用、思考与练习38
6.4.1.计算最佳水槽断面面积38
6.4.2对约束优化的讨论39
6.4.3.工程优化问题的计算39
实验7.隐函数、方程求根、不动点和迭代39
7.1知识要点与背景39
7.1.1隐函数存在定理与四连杆机构的运动40
7.1.2不动点和函数迭代40
7.2实验与观察40
7.2.1隐函数的存在定理的可视化40
7.2.2.用蛛网图观察不动点迭代41
7.2.3简单和复杂：
二次函数的迭代和混沌43
7.3.应用、思考与练习44
7.3.1四连杆输出角的运动规律和动画模拟44
7.3.2轨道飞行器的拦截45
7.3.3怎样保证或加速迭代序列的收敛45
7.3.4.混沌有哪些特点？
46
7.4非线性方程组求解46
实验8.线性代数实验46
8.1实验（Ⅰ）：
用Matlab学线性代数46
8.1.1实验与观察：
向量组的线性关系和解线性方程组47
8.1.2应用、思考与练习49
8.2实验（Ⅱ）：
矩阵的相似化简50
8.2.1实验与观察：
矩阵的特征―相似标准形的作用50
8.2.2应用、思考与练习51
实验9概率统计实验53
9.1实验（I）：
Galton钉板试验53
9.1.1实验与观察:
Galton钉板模型和二项分布53
9.1.2应用、思考和练习55
9.2实验（Ⅱ）：
热轧机的调整56
9.2.1实验与观察:
控制粗轧的浪费57
9.2.2应用、思考与练习59
9.3实验（Ⅲ）参数估计和假设检验60
9.3.1实验与观察:
极大似然估计61
9.3.2应用、思考与练习62
实验10.数值仿真66
10.1知识要点与背景：
单自由度阻尼系统66
10.2.2振动弹簧的实时动画67
10.3.3物理问题的数值模拟68
实验11.傅氏分析与小波分析69
11.1知识要点—傅氏分析与小波分析69
11.1.1傅氏分析69
11.1.2小波分析69
11.2观察与实验69
11.2.1信号频谱分析69
11.2.2如何得到小波函数70
11.2.3单尺度一维离散小波分解与重构71
11.2.4多尺度分解与重构71
11.3应用、思考与练习72
11.3.1信号的奇异性检测72
11.3.2信号去噪73
11.3.3信号的压缩73
11.3.4练习73
实验12．金融分析实验73
12.1知识要点和背景：
最优投资组合及其计算73
12.1.1包含无风险证券的投资组合73
12.1.2无风险证券投资组合的计算73
12.2机会的价值73
12.2.1简单二项式模型机会价值74
12.2.2两期二项式模型机会价值74
12.2.3观察与思考75
《数学实验与Matlab》程序
周晓阳
华中科技大学数学系
我将程序按书中的顺序排列，这样便于读者利用。
本书程序均通过了调式。
可直接拷贝到命令窗口运行或编制M文件运行。
如出现问题，可能是中英文标点的缘故（排版有可能使用了中文标点），请将中文标点换为英文标点试试。
实验１．矩阵运算与Matlab命令
1.1知识要点与背景：
知识要点和背景：
日常矩阵及其运算
【A=[423;132;133;322],%表1-1、表1-2的数据分别写成矩阵形式
B=[35206045;10155040;20124520]】
【C=A*B%矩阵乘法，求各订单所对应的原材料和劳动力。
】
【whos%查看Matlab工作空间中变量及其规模】
1.2实验与观察：
矩阵和Matlab语言
1.2.1向量的生成和运算
【x=linspace（0,4*pi,100）;%将[0,4π]区间100等分，产生了一个100维向量
y=sin（x）;%计算函数值，产生了一个与x同维的100维函数向量y
y1=sin（x）.^2;%计算函数向量，注意元素群运算
y2=exp（-x）.*sin（x）;
%以x为横坐标，y为纵坐标画函数的图用不同的线型将函数曲线绘制在一个图上
plot（x,y,'-',x,y1,'-',x,y2,'.-'）】
1.向量的创建
◆直接输入向量。
【x1=[124],x2=[1,2,1],x3=x1'】
◆冒号创建向量。
【x1=3.4:
6.7
x2=3.4:
2:
6.7
x3=2.6:
-0.8:
0】
◆生成线性等分向量。
【x=linspace（0,1,5）】
2.向量的运算
【y=sin（x）】
【y1=sin（x）.^2;y2=exp（-x）.*sin（x）;】
1.2.2.矩阵创建和运算
1.创建矩阵
（1）数值矩阵的创建
◆直接输入法创建简单矩阵。
【A=[1234;5678;9101112]】
【B=[-1.3,sqrt（3）;（1+2）*4/5,sin（5）;exp
（2）,6]】
（2）符号矩阵的创建
◆
【symsa11a12a13a14a21a22a23a24a31a32a33a34…
b11b12b13b14b21b22b23b24b31b32b33b34
A1=[a11a12a13a14;a21a22a23a24;a31a32a33a34],
B1=[b11b12b13b14;b21b22b23b24;b31b32b33b34]】
2.矩阵的运算
【C=A1+B1,D=A1-B1】
【symsc
cA=c*A1】
【C=A1*B1】
{?
?
?
Errorusing==>sym/mtimes,Innermatrixdimensionsmustagree.}
【A2=A1（:
1:
3）,B1】
【G=A2*B1】
【g11=A2（1,:
）*B1（:
1）】
【A,A_trans=A'】
【H=[123;210;123],K=[123;210;231]
h_det=det（H）,k_det=det（K）,H_inv=inv（H）,K_inv=K^-1】
【A=[301;110;014];
B=inv（A-2*eye（3））*A,B=（A-2*eye（3））\A】
3.分块矩阵：
矩阵的裁剪、分割、修改与抽取
（1）
【A=[10112;01-123;30510;23121],vr=[1,3];vc=[1,3];
A1=A（vr,vc）%取出A的1、3行和1、3列的交叉处元素构成新矩阵A1】
◆将上面的矩阵A分为四块，并把它们赋值到矩阵B中，观察运行后的结果。
【A11=A（1:
2,1:
2）,A12=A（1:
2,3:
5）,A21=A（3:
4,1:
2）,A22=A（3:
4,3:
5）
B=[A11A12;A21A22]】
A=
205
421
0-12
B=
124-1
5310
-1023
C=
-341813
131420-1
-7-336
（2）矩阵的修改和提取
◆【A=[10112;01-123;30510;23121]
A（1,:
）=[00000];A】
◆观察：
【B（:
[2,4]）=[]%删除矩阵B的第2、4列】
（3）矩阵元素的抽取
4.生成特殊矩阵
。
◆
【y1=rand（1,5）,y2=rand（1,5）,
rand（'seed',3）,x1=rand（1,5）,rand（'seed',3）,x2=rand（1,5）】
5.常用矩阵函数
6.数据的简单分析
◆
【rand（'seed',1）;A=rand（3,6）,
Asort=sort（A）,Amax=max（A）,Asum=sum（A）】
1.2.3Matlab工作环境和编程
2.Matlab的基本设计
1.3应用、思考与练习
1.3.1关系矩阵
1.3.2投入产出
1.3.3循环比赛的名次
【A=[0110;0011;0001;1000],
e=ones（4,1）;c=A*e;
s=c'】
★画矩阵结构图的gplot指令。
◆（3）
【clf,A=[0110;0011;0001;1000];xy=[01;00;-1–0.5;1–0.5];
graphy_plot（A,xy,1,0.5）,%gplot（A,xy）】
1.3.4参考程序
graphy_plot.m
【functiony=graphy_plot（A,xy,l,p）
%画矩阵的有向结构图。
A为邻接矩阵，xy为顶点坐标，l控制参数，l=0，画无向图；
%l~=0,画有向图。
p为控制箭头大小的参数。
a=-max（abs（xy（:
1）））*1.1;b=max（abs（xy（:
1）））*1.1;
c=-max（abs（xy（:
2）））*1.1;d=max（abs（xy（:
2）））*1.1;
ifl=0
gplot（A,xy）,axis（[abcd]）,holdon,
elseifl~=0
U=[];V=[];X=[];Y=[];
n=length（A（:
1））;
fori=1:
n
k=find（A（i,:
）~=0）;
m=length（k）;
if（m~=0）
forj=1:
m
u
（1）=（xy（k（j）,1）-xy（i,1））;v
（1）=（xy（k（j）,2）-xy（i,2））;
u
（2）=eps;v
（2）=eps;U=[u;U];V=[v;V];
X=[[xy（i,1）xy（k（j）,1）];X];Y=[[xy（i,2）xy（k（j）,2）];Y];
end
text（xy（i,1）,xy（i,2）,['\bullet\leftarrow\fontsize{16}\it{V}',…
um2str（i）]）;holdon,
end
end
gplot（A,xy）,axis（[abcd]）,holdon,
h=quiver（X,Y,U,V,p）;set（h,'color','red'）;holdon,
plot（xy（:
1）,xy（:
2）,'k.','markersize',12）,holdon,
end,holdoff】
实验２．函数的可视化与Matlab作
2.1实验与观察：
函数的可视化
2.1.1Matlab二维绘图命令
1.周期函数与线性p－周期函数
◆观察：
【clf,x=linspace（0,8*pi,100）;
F=inline（'sin（x+cos（x+sin（x）））'）;
y1=sin（x+cos（x+sin（x）））;y2=0.2*x+sin（x+cos（x+sin（x）））;
plot（x,y1,'k:
',x,y2,'k-'）legend（'sin（x+cos（x+sin（x））','0.2x+sin（x+cos（x+sin（x）））',2）】
2.plot指令：
绘制直角坐标的二维曲线
3.图形的属性设置和屏幕控制
【h=plot（[0:
0.1:
2*pi],sin（[0:
0.1:
2*pi]））;gridon
set（h,'LineWidth',5,'color','red'）;set（gca,'GridLineStyle','-','fontsize',16）】
◆设置y坐标的刻度并加以说明，并改变字体的大小。
【h=plot（[0:
0.1:
2*pi],sin（[0:
0.1:
2*pi]））;gridon，
set（gca,'ytick',[-1-0.500.51]）,set（gca,'yticklabel','a|b|c|d|e'）,
set（gca,'fontsize',20）】
4.文字标注指令
【plot（x,y1,'b',x,y2,'k-'）,
set（gca,'fontsize',15,'fontname','timesNewRoman'）,%设置轴对象的字体为times
%NewRoman，字体的大小为15
title（'\it{Peroidandlinearperoidfunction}'）;%加标题，注意文字用单引号''加上
%斜杠'\'后可输入不同的设置，例如it{…}表示花括号里的文字为斜体；如果有多项设置，
%则可用\…\…\…连续输入。
xlabel（'xfrom0to8*piit{t}\'）;ylabel（'\it{y}'）;%说明坐标轴
text（x（49）,y1（50）-0.4,'\fontsize{15}\bullet\leftarrowTheperiodfunction{\itf（x）}'）;
%在坐标（x（49）,y1（50）-0.4）处作文字说明，各项设置用"\"隔开。
%\fontsize{15}\bullet\leftarrow的意义依次是：
\字体大小=15\画圆点\左箭头
text（x（14）,y2（50）+1,'\fontsize{15}Thelinearperiodfunction{\itg（x）}\rightarrow
\bullet'）%与上一语句类似，用右箭头】
图2.5文字标注
◆观察指令legend和num2str的用法：
在同一张图上画出
，这里
并进行适当的标注。
zxy2_2.m
【clf,t=0:
0.1:
3*pi;alpha=0:
0.1:
3*pi;
plot（t,sin（t）,'r-'）;holdon;plot（alpha,3*exp（-0.5*alpha）,'k:
'）;
set（gca,'fontsize',15,'fontname','timesNewRoman'）,
xlabel（'\it{t（deg）}'）;ylabel（'\it{magnitude}'）;
title（'\it{sinewaveand{\it{Ae}}^{-\alpha{\itt}}wave}'）;%注意\alpha的意义
text（6,sin（6）,'\fontsize{15}TheValue\it{sin（t）}at{\itt}=6\rightarrow\bullet','HorizontalAlignment','right'）,
%上面的语句是一整行，如果要写成两行，必须使用续行号…，例如要在“bullet',”
%后换行，需写“bullet',…”后才能换行。
%'HorizontalAlignment','right'表示箭头所指的曲线对象在文字的右边。
text（2,3*exp（-0.5*2）,['\fontsize{15}\bullet\leftarrowTheValueof\it{3e}^{-0.5\it{t}}=',num2str（3*exp（-0.5*2））,'at\it{t}=2']）;
%num2str的用法：
['string1',num2str,'string2']，注意方括号的使用。
%legend（'\itsin（t）','{\itAe}^{-\alphat}'）%请结合图形观察此命令的使用】
运行结果如图2.6所示。
5.图形窗口的创建和分割
◆观察：
【clf，b=2*pi;x=linspace（0,b,50）;
fork=1:
9
y=sin（k*x）;
subplot（3,3,k）,plot（x,y）,axis（[0,2*pi,-1,1]）
end】
x=linspace（0,4*pi,150）;
y1=cos（x）;
y=cos（cos（cos（cos（cos（cos（cos（x）+sin（x）））+sin（x）））+sin（x）））;
plot（x,y,x,y1,'r-'）,grid
2.1.2多元函数的可视化与空间解析几何（三维图形）
本节通过高等数学的几个例子观察Matlab的三维绘图功能和技巧。
1.绘制二元函数
◆观察：
绘制
的图象，作定义域的裁剪。
◆
（1）观察meshgrid指令的效果。
【a=-0.98;b=0.98;c=-1;d=1;n=10;
x=linspace（a,b,n）;y=linspace（c,d,n）;
[X,Y]=meshgrid（x,y）;
plot（X,Y,'+'）】
★三维绘图指令mesh、meshc、surf。
◆
（2）做函数的定义域裁剪，观察上述三维绘图指令的效果。
程序zxy2_4.m
【clear,clf,
a=-1;b=1;c=-15;d=15;n=20;eps1=0.01;
x=linspace（a,b,n）;y=linspace（c,d,n）;
[X,Y]=meshgrid（x,y）;
fori=1:
n%计算函数值z，并作定义域裁剪
forj=1:
n
if（1-X（i,j））z（i,j）=NaN;%作定义域裁剪，定义域以外的函数值为NaN
else
z（i,j）=1000*sqrt（1-X（i,j））^-1.*log（X（i,j）-Y（i,j））;
end
end
end
zz=-20*ones（1,n）;plot3（x,x,zz）,gridoff,holdon%画定义域的边界线
mesh（X,Y,z）%绘图，读者可用meshz,surf,meshc在此替换之
xlabel（'x'）,ylabel（'y'）,zlabel（'z'）,boxon%把三维图形封闭在箱体里】
◆运行了zxy2_4.m以后，有关向量存储在工作空间中，在此基础上，观察上述等值线绘制指令的运行效果。
【[cs,h]=contour（X,Y,z,15）;clabel（cs,h,'labelspacing',244）】
2.三元函数可视化:
slice指令
◆观察：
绘制三元函数
的可视化图形。
【clf,x=linspace（-2,2,40）;y=x;z=x;
[X,Y,Z]=meshgrid（x,y,z）;w=X.^2+Y.^2+Z.^2;
slice（X,Y,Z,w,[-1,0,1],[-1,0,1],[-1,0,1]）,colorbar】
3.空间曲线及其运动方向的表现：
plot3和quiver指令
【clf,t=0:
0.1:
1.5;
Vx=2*t;Vy=2*t.^2;Vz=6*t.^3-t.^2;
x=t.^2;y=（2/3）*t.^3;z=（6/4）*t.^4-（1/3）*t.^3;%由速度得到曲线
plot3（x,y,z,'r.-'）,holdon%画飞行轨迹
%算数值梯度，也就是重新计算数值速度矢量，这只是为了编程的方便，不是必须的
图2.12飞机的飞行轨迹与方向
Vx=gr
升级会员 