MATLAB数学实验全套课件电子教案板

MATLAB数学实验2024/9/152第一章Matlab入门第一章MATLAB入门1.1MATLAB桌面1.2数据和变量1.3数组及其运算1.4字符串、元胞和结构

2024/9/153第一章Matlab入门1.1MATLAB桌面最小安装：MATLAB2013a，CurveFittingToolbox，OptimizationToolbox，SymbolicMathToolbox工具条：Home,Plots,Apps等窗口:指令(CommandWindow)当前文件夹(CurrentDirectory)工作空间(Workspace)指令历史(CommandHistory)例:a=1;b=2;c=a+b(不输入提示符>>)2024/9/154第一章Matlab入门1.2数据和变量P4例1.1球的体积计算表达式分号(;),逗号(,),省略号(...)。历史指令调用科学记数法数据显示格式Shortlongrational显示格式与计数精度区别2024/9/155第一章Matlab入门1.2数据和变量复数i,j注意(-8)^(1/3)预定义变量pi圆周率3.1415…

eps浮点数识别精度2.22×10-16realmin最小正实数2.2251×10-308realmax最大正实数1.7977×10308

Inf无穷大NaN不定值注意e在MATLAB中只是普通变量

2024/9/156第一章Matlab入门1.2数据和变量用户变量命名规则：字母开头，由字母、数字或下划线组成（不得使用减号、空格等

），区分大小写防止与系统的预定义变量名(如i,j,pi,eps等)，函数名(如who,length等)，保留字(for,if,while,end等)冲突。特殊变量ans

是系统本身一个特殊变量名，若运算结果没有赋于任何变量，系统将其赋予ansclear

清除(注意ClearWorkspace与ClearCommandWindow的区别.)2024/9/157第一章Matlab入门1.2数据和变量数据Mat文件实现与外部数据文件交换:mat,txt等菜单方式：SaveWorkspace和ImportData例:save-clear-import指令方式：save和loadMATLAB还提供了方便的工具来导入外部数据文件，包括文本文件、Excel(Spreadsheets)文件、图像文件等，以便与其它应用程序交换数据，我们将在第二章介绍。

2024/9/158第一章Matlab入门1.3数组及其运算a=[123;456;780]变量编辑器(Variables)和剪贴板的使用数组的输入和分析中括号[]表示矩阵，同行无素间用空格或逗号分隔，不同行间用分号或回车分隔。冒号运算函数linspace(x1,x2,n)生成x1与x2间的n维等距行向量,即将[x1,x2]n-1等分length,size编址：不能为0，按列编址,如a(6)2024/9/159第一章Matlab入门1.3数组及其运算数组的抽取和统计

分块矩阵a([1,3],1:3)sum,prod,min,max按列计算例如[x,y]=max([-13;5-6;24])x=54y=232024/9/1510第一章Matlab入门1.3数组及其运算数组运算

A+B与A-B加与减

k*A或A*k数乘矩阵

k+A与k-Ak加(减)A的每个元素A.^k,k.^A数组乘方

A.*B数组乘数组

k./A数除以数组左除A.\B=右除B/.A数组除法

点运算就是对应元素的运算P12例（注意点运算与矩阵运算的区别）2024/9/1511第一章Matlab入门1.3数组及其运算数学函数矩阵的数学函数也是按元素的运算，使用通常的函数号，如sin(A),cos(A),asin(A),acos(A),tan(A),cot(A),exp(A),sqrt(A)等。fix向0取整floor向下取整ceil向+取整mod模余rem除法余数abs绝对值(模）real复数实部imag复数虚部angle复数幅角conj复数共轭log

自然对数ln

log10以10为底对数2024/9/1512第一章Matlab入门1.3数组及其运算4、关系与逻辑运算

<、<=小于、小于等于

>、>=大于、大于等于

==、~=等于、不等于&（与）、|（或）、~（非）any、all、find“真(true)”用1表示，“假(false)”用0.逻辑运算中，所有非零元素作为1处理.P15例子小结：Matlab中的三类运算函数按元素运算（数学函数）：如sin(A),log(A),find(A)按列运算（统计函数）：如sum(A),prod(A),min(A),max(A),any(A),all(A)按矩阵运算（矩阵函数）：如size(A),length(A),det(A)2024/9/1514第一章Matlab入门1.4字符串、元胞和结构数据类型：数值(Double)逻辑(Logical)字符(Char)

元胞(Cell)

结构(Structure)2024/9/1515第一章Matlab入门1.4字符串、元胞和结构字符串（P17例）单引号(英文半角输入状态！)中文字符不要在word中输入后copy引号内字符显示应为淡紫色字符串拼接字符串转化double,char,num2str,str2num比较：a=’12’,b=double(a),c=str2num(a)Eval执行字符串书写的指令(P18例)2024/9/1516第一章Matlab入门1.4字符串、元胞和结构

元胞和结构（P18例）字符矩阵char数值与字符混合元胞{}结构:域的概念struct2cell和cell2struct2024/9/1517第一章Matlab入门习题P20ex1,ex2,ex3,ex4MATLAB数学实验第二章MATLAB编程与作图2024/9/1519第二章MATLAB编程与作图第二章MATLAB编程与作图2.1程序设计2.2作图

2.3在线帮助和文件管理2024/9/1520第二章MATLAB编程与作图2.1程序设计循环语句for循环变量=初值:增量:终值,语句；endwhile（条件式）,语句；end分支语句if（条件式）,语句;endif（条件式1),语句1;elseif（条件式2),

语句2;……;else,语句;endswitch(分支变量)case(值1),语句1;case(值2),语句2;……;otherwise语句;end其它：pause,break,return,error2024/9/1521第二章MATLAB编程与作图2.1程序设计»s=0;forn=1:100,s=s+1/n/n;end；s»clear;s=0;n=1;whilen<=100,s=s+1/n/n;n=n+1;end;s强行中断：Ctrl+C

2024/9/1522第二章MATLAB编程与作图2.1程序设计M脚本文件eg2_1在Editor窗口文件名一律以字母开头，以字母、数字或下划线组成。注意不要含有空格、减号等.!!!常见错误：用1或eg2-1作为文件名

M文件名一般都用小写字母注意：必须保存在当前目录(CurrentDirectory)或matlab路径中2024/9/1523第二章MATLAB编程与作图2.1程序设计M函数文件function输出变量=函数名(输入变量)

语句;eg2_1fM函数调用必须给予输入参数值M函数中变量为局部变量2024/9/1524第二章MATLAB编程与作图2.1程序设计M函数文件M函数在edit窗口编写，在command窗口调用!!!常见错误：将function写在commond窗口M函数是以该函数的磁盘文件主名调用，而不是文件中的函数名称!!!常见错误：错误用文件中的函数名称调用，却与磁盘文件主名不同2024/9/1525第二章MATLAB编程与作图2.1程序设计函数句柄(handle)编写M函数eg2_1f.m函数句柄fname=@eg2_1f;feval(fname,1000)或者fname(1000)注意feval与eval的区别>>m=1000;>>str='sum(1./(1:m).^2';>>eval(str)2024/9/1526第二章MATLAB编程与作图inline函数与匿名函数Inline函数(即将被淘汰)fun=inline(‘函数表达式’,‘自变量’)>>fname=inline('sum(1./(1:m).^2)','m')>>fname(1000)匿名函数fun=@(自变量)函数表达式匿名函数的参数传递更灵活>>k=2;fname=@(m)sum(1./(1:m).^k)>>fname(1000)2024/9/1527第二章MATLAB编程与作图2.1程序设计注释：%开头，对本行后面字符起作用，不参与运算。

对话：input，disp

全程变量与局部变量

nargin、nargout和varargin

子函数和嵌套函数提高速度2024/9/1528第二章MATLAB编程与作图2.1程序设计普通编程functions=f(m)s=0;forn=1:ms=s+1/n/n;end向量化编程functions=f(m)n=1:m;s=sum(1./n.^2);尽量少用for语句2024/9/1529第二章MATLAB编程与作图本书配套工具箱的使用主要内容教材所有程序；数学实验常用Matlab指令的中文帮助安装步骤1.将ME.rar解压缩至matlab\toolbox\;2.启动Matlab,利用File菜单中的Setpath将matlab\toolbox\me增至path中,放在最底部，并保存设置；3.回到你的工作目录。现在ME已成为一个普通的工具箱了。2024/9/1530第二章MATLAB编程与作图2.1程序设计

例2.4编一M函数,对任意输入的向量x,可计算分段函数值构成的向量。分量方式eg2_4a，慢向量方式eg2_4b,eg2_4c,快数组预分配y=zeros(size(x))2024/9/1531第二章MATLAB编程与作图2.2作图曲线图plot(x,y)

以数据(x(i),y(i))为节点的折线图,其中x,y为同长度的向量

plot(x1,y1,x2,y2,...)

多组数据折线图

fplot(fun,[a,b])

函数fun在区间[a,b]上的函数图

plot3(x,y,z)

空间曲线图,其中x,y,z为同长度的向量图形导出到word线型与标记P31表eg2_5

曲线图y=x3-x-1和y=|x|0.2sin(5x)

2024/9/1532第二章MATLAB编程与作图2.2作图曲面图[x,y]=meshgrid(xa,ya)当xa,ya分别为m维和n维行向量，得到x和y均为n行m列矩阵。meshgrid常用于生成X-Y平面上的网格数据。[x,y]=ndgrid(xa,ya)与meshgrid类似，但得到的x和y均为m行n列矩阵

mesh(x,y,z)

绘制网面图，是最基本的曲面图形命令,其中x,y,z是同阶矩阵，表示曲面三维数据。surf(x,y,z)

绘制曲面图，与mesh用法类似。2024/9/1533第二章MATLAB编程与作图meshgrid解释xa=6:8;ya=1:4;[x,y]=meshgrid(xa,ya)%生成X-Y面上网格z=x.^2+y.^2%计算XY面各格点的z轴高度xyz6781113750656782224053686783334558736784445265802024/9/1534第二章MATLAB编程与作图ngrid解释xa=6:8;ya=1:4;[x,y]=ngrid(xa,ya)%生成X-Y面上网格z=x.^2+y.^2%计算XY面各格点的z轴高度xyz666612343740455277771234505358658888123465687380eg2_6二元函数图z=xexp(-x2-y2)2024/9/1535第二章MATLAB编程与作图2024/9/1536第二章MATLAB编程与作图2.2作图图形说明和定制title标题说明；xlabel，ylabel，zlabel说明坐标轴x,y,z;holdon/holdoff

保留/释放现有图形axis([a,b,c,d])确定坐标轴范围a<x<b,c<y<daxis([a,b,c,d,e,f])定制3维坐标轴范围figure\close开\关一个新图形窗口subplot(m,n,k)

将图形窗口分为m*n个子图,指向第k幅图legend(str1,str2,...)图例eg2_7空间曲线2024/9/1537第二章MATLAB编程与作图2.2作图图形窗口菜单和工具栏

图形编辑2024/9/1538第二章MATLAB编程与作图2.3在线帮助和文件管理在线帮助helphelp子目录名help命令或函数doc命令或函数(超文本帮助)lookfor关键字docsearch关键字

(超文本搜索)

typeM文件主名whichM文件主名2024/9/1539第二章MATLAB编程与作图2.3在线帮助和文件管理文件和文件夹管理MATLAB接受到一个命令的搜索过程初学者在M文件的保存上常出现几种错误

设置你自己的工作目录

(Currentdirectory)

设置MATLAB默认搜索路径(Path)队列2024/9/1540第二章MATLAB编程与作图2.3在线帮助和文件管理读写外部数据文件save存为Mat数据文件loadMat数据文件载入xlswrite数据写入Excel文件xlsread读取Excel文件dlmwrite,fprintf写入文本文件dlmread,textscan,textread,fscanf读取文本文件imread读取图像文件importdata导入各种数据文件工具条Home:Importdata2024/9/1541第二章MATLAB编程与作图2.3在线帮助和文件管理P42例：Excel文件的调用>>[Grade,Name,Raw]=xlsread('c:\grade.xls','Sheet1','a2:d3')>>Total=Grade(:,1)+Grade(:,2)+Grade(:,3);>>xlswrite('c:\grade.xls',Total,'Sheet1','e2:e3')>>xlswrite('c:\grade.xls',{'总分'},'Sheet1','e1')2024/9/1542第二章MATLAB编程与作图习题P46ex2,ex3,ex5,ex6(单号)MATLAB数学实验第三章

矩阵代数

2024/9/1544第三章矩阵代数第三章

矩阵代数3.1预备知识：线性代数3.2矩阵代数的MATLAB指令3.3计算实验：线性方程组的通解3.4建模实验：投入产出分析和基因遗传2024/9/1545第三章矩阵代数3.1预备知识：线性代数线性方程组记为Ax=b2024/9/1546第三章矩阵代数3.1预备知识：线性代数线性方程组若秩(A)

秩(A,b)，则无解；若秩(A)=秩(A,b)=n,存在唯一解；若秩(A)=秩(A,b)<n,存在无穷多解；通解是齐次线性方程组Ax=0的基础解系与Ax=b的一个特解之和。2024/9/1547第三章矩阵代数3.1预备知识：线性代数逆矩阵方阵A称为可逆的，如果存在方阵B，使AB=BA=E，记B=A-1方阵A可逆的充分必要条件:

A

0A-1=A*/|A|这里A*为A的伴随矩阵（AE）行变换（EA-1）2024/9/1548第三章矩阵代数3.1预备知识：线性代数特征值与特征向量对于方阵A，若存在数

和非零向量x使Ax=

x，则称

为A的一个特征值，x为A的一个对应于特征值

的特征向量。特征值计算归结为特征多项式的求根。特征向量计算：齐次线性方程组

(A-

E)x=0的所有一组线性无关解。2024/9/1549第三章矩阵代数3.2矩阵代数的MATLAB指令运算符A’(共轭)转置,A.’

转置A+B与A-B

加与减k+A与k-A

数与矩阵加减k*A或A*k

数乘矩阵 A*B

矩阵乘法A^k

矩阵乘方左除A\B

为AX=B的解右除B/A

为XA=B的解 与数组运算不同与数组运算相同2024/9/1550第三章矩阵代数3.2矩阵代数的MATLAB指令矩阵运算与数组运算的区别数组运算按元素定义，矩阵运算按线性代数定义矩阵的加、减、数乘等运算与数组运算是一致的

注意：矩阵的乘法、乘方、除法与数组乘法、乘方、除法不同！数与矩阵加减、矩阵除法在数学上是没有意义的。但在MATLAB中有定义。

例子P51-522024/9/1551第三章矩阵代数3.2矩阵代数的MATLAB指令特殊矩阵生成zeros(m,n)m行n列的零矩阵;ones(m,n)m行n列的元素全为1的阵;eye(n)n阶单位矩阵;rand(m,n)m行n列[0,1]上均匀分布随机数矩阵举例说明2024/9/1552第三章矩阵代数3.2矩阵代数的MATLAB指令矩阵处理

trace(A)迹(对角线元素的和)diag(A)

A对角线元素构成的向量;diag(x)

向量x的元素构成的对角矩阵.tril(A)A的下三角部分triu(A)A的上三角部分flipud(A)矩阵上下翻转fliplr(A)矩阵左右翻转reshape(A,m,n)

矩阵A的元素重排成m行n列矩阵

2024/9/1553第三章矩阵代数3.2矩阵代数的MATLAB指令矩阵分析

rank(A)

秩det(A)

行列式;inv(A)

逆矩阵;null(A)

Ax=0的基础解系；null(A,’r’)基础解系的有理(整数或分数)形式orth(A)

A列向量正交规范化norm(x)向量x的范数（长度，模）norm(A)矩阵A的范数2024/9/1554第三章矩阵代数3.2矩阵代数的MATLAB指令特征值与标准形p=poly(A)

返回方阵A的特征多项式eig(A)

方阵A的特征值[V,D]=eig(A)返回方阵A的特征值和特征向量。其中D为的特征值构成的对角阵，每个特征值对应的V的列为属于该特征值的一个特征向量。[V,J]=jordan(A)返回A的相似变换矩阵和约当标准形(A的相似对角化)

例子P56-572024/9/1555第三章矩阵代数3.3计算实验：线性方程组求解

矩阵除法

(1)当A为方阵，A\B结果与inv(A)*B一致；(2)当A不是方阵,AX=B存在唯一解,A\B将给出这个解；(3)当A不是方阵,AX=B为不定方程组(即无穷多解)，A\B将给出一个具有最多零元素的特解；(4)当A不是方阵,AX=B若为超定方程组（即无解）,A\B给出最小二乘意义上的近似解，即使得向量AX－B的范数达到最小。

2024/9/1556第三章矩阵代数3.3计算实验：线性方程组求解例3.1解方程组

2024/9/1557第三章矩阵代数3.3计算实验：线性方程组求解例3.2线性方程组通解用rref化为行最简形以后求解用除法求出一个特解，再用null求得一个齐次组的基础解系用符号数学工具箱中的solve求解(第七章)

2024/9/1558第三章矩阵代数3.3计算实验：线性方程组求解相似对角化及应用

如果n阶方阵A有n个线性无关的特征向量，则必存在正交矩阵P,使得P-1AP=

,其中

是A的特征值构成的对角矩阵，P的列向量是对应的n个正交特征向量。使用MATLAB函数eig求得的每个特征向量都是单位向量(即模等于1)，并且属于同一特征值的线性无关特征向量已正交化，所以由此容易进行相似对角化。

2024/9/1559第三章矩阵代数3.3计算实验：线性方程组求解例3.3

用相似变换矩阵P将A相似对角化，并求>>A=[11/40;01/20;01/41];[P,T]=eig(A)AP=PT

P-1AP=T

A=PTP-1An=PTP-1PTP-1……PTP-1=PTnP-12024/9/1560第三章矩阵代数3.4建模实验设有n个经济部门，xi为部门i的总产出，cij为部门j单位产品对部门i产品的消耗，di为外部对部门i的需求，fj为部门j新创造的价值。分配平衡方程组(部门i产品=内部需求+外部需求)消耗平衡方程组(部门j产值=生产消耗+新创造价值)

2024/9/1561第三章矩阵代数列昂杰夫Leontief,Nobel经济学奖

1973年W.W.LeontiefwasaRussian-Americaneconomistnotableforhisresearchonhowchangesinoneeconomicsectormayhaveaneffectonothersectors.LeontiefwontheNobelMemorialPrizeinEconomicSciencesin1973,andthreeofhisdoctoralstudentshavealsobeenawardedtheprize(PaulSamuelson1970,RobertSolow1987,VernonSmith2002).2024/9/1562第三章矩阵代数投入产出分析令C=（cij），X=(x1,…,xn)’,D=(d1,…,dn)’，F=(f1,…,fn)’,则分配平衡方程组X=CX+D令A=E－C，E为单位矩阵，则

AX=DC称为直接消耗矩阵A称为列昂杰夫矩阵。2024/9/1563第三章矩阵代数Y=[1,1,…,1]B（B各列的和）Y表示各部门的总投入(消耗)。新创造价值F=X–Y'B=CB表示各部门间的投入产出关系，称为投入产出矩阵。注：bij=cijxj2024/9/1564第三章矩阵代数投入产出关系表(行：分配平衡，列：消耗平衡)

消耗部门外界需求总产出123生产部门1b11b12b13d1x12b21b22b23d2x23b31b32b33d3x3新创造价值f1f2f3

总产出x1x2x3注：bij=cijxj2024/9/1565第三章矩阵代数投入产出平衡行：分配平衡X=B各行之和+外界需求列：消耗平衡X=B各列之和+新创造价值2024/9/1566第三章矩阵代数投入产出分析

例3.4某地有三个产业，一个煤矿，一个发电厂和一条铁路，开采一元钱的煤，煤矿要支付0.25元的电费及0.25元的运输费;

生产一元钱的电力，发电厂要支付0.65元的煤费，0.05元的电费及0.05元的运输费;

创收一元钱的运输费,铁路要支付0.55元的煤费和0.10元的电费，在某一周内煤矿接到外地金额50000元定货，发电厂接到外地金额25000元定货，外界对地方铁路没有需求。2024/9/1567第三章矩阵代数解：这是一个投入产出分析问题。设x1为本周内煤矿总产值，x2为电厂总产值,x3为铁路总产值,则问三个企业间一周内总产值多少才能满足自身及外界需求？三个企业间相互支付多少金额？三个企业各创造多少新价值?2024/9/1568第三章矩阵代数直接消耗矩阵C=外界需求向量D=产出向量X=则原方程为(E-C)X=D投入产出矩阵为

B=C*diag(X)总投入向量

Y=ones(1,3)*B新创造价值向量

F=X-Y’2024/9/1569第三章矩阵代数表3.3投入产出分析表(单位：元)

消耗部门外界需求总产出煤矿电厂铁路生产部门煤矿0365061558250000102088电厂25522280828332500056163铁路2552228080028330新创造价值51044140419915

总产出10208856163283302024/9/1570第三章矩阵代数例5设金鱼某种遗传病染色体的正常基因为A，异常基因为a,那么AA，Aa，aa分别表示正常金鱼，隐性患者，显性患者。设初始分布为90%正常金鱼，10%的隐性患者，无显性患者。考虑下列两种配种方案对后代该遗传病基因型分布的影响方案一：同类基因结合，均可繁殖；方案二：显性患者不允许繁殖，隐性患者必须与正常金鱼结合繁殖2024/9/1571第三章矩阵代数后代基因型的概率后代是从父母体的基因对中各继承一个基因，形成自己的基因型。2024/9/1572第三章矩阵代数第一种方案同类基因结合，均可繁殖；2024/9/1573第三章矩阵代数第二种方案显性患者不允许繁殖，隐性患者必须与正常金鱼结合繁殖2024/9/1574第三章矩阵代数基因遗传矩阵表示2024/9/1575第三章矩阵代数解设初始分布X(1)=(0.90.10)’,第n代分布为X(n)=方案一：X(n)=Mn-1X(1)方案二：X(n)=Nn-1X(1)2024/9/1576第三章矩阵代数结论计算方法：取n足够大，如n=20方案一：Mn-1X(1)

[0.95,0,0.05]’存在5%的显性患者方案二：Nn-1X(1)

[1,0,0]’显性患者和隐性患者都消失本例说明了杂交的优势另一计算方法：用矩阵相似对角化2024/9/1577第三章矩阵代数习题P65ex2,ex3,ex4,ex5,ex6,ex102024/9/1578第三章矩阵代数习题4转移矩阵表示2024/9/1579第三章矩阵代数习题5投入产出关系

消耗部门外界需求总产出123生产部门1c11x1c12x2c13x3d1x12c21x1c22x2c23x3d2x23c31x1c32x2c33x3d3x3

注：bij=cijxjMATLAB数学实验第四章函数和方程第四章函数和方程4.1预备知识：零点、极值和最小二乘法4.2函数零点、极值和最小二乘拟合的MATLAB指令4.3计算实验：迭代法4.4建模实验：购房贷款的利率4.1预备知识：零点非线性方程

f(x)=0若对于数

有f(

)=0,则称

为方程的解或根，也称为函数f(x)的零点非线性方程求解通常用数值方法求近似解.非线性方程（组）f(x)=0，其中x=(x1,x2,…,xn),f=(f1,f2,…,fm)

4.1预备知识：极值设x为标量或向量，y=f(x)是x

D上的标量值函数。如果对于包含x=a的某个邻域

,有f(a)f(x)

(f(a)f(x))对任意x

成立,则称a为f(x)的一个局部极小(大)值点。。如果对任意x

D，有f(a)f(x)（f(a)f(x)）成立，则称a为f(x)在区域D上的一个全局极小(大)值点。4.1预备知识：极值4.1预备知识：最小二乘拟合4.1预备知识：最小二乘拟合成绩选手国籍日期13秒24米尔布恩美国1972年9月7日13秒21卡萨那斯古巴1977年8月21日13秒16内赫米亚美国1979年4月14日13秒00内赫米亚美国1979年5月6日12秒93内赫米亚美国1981年8月19日12秒92金多姆美国1989年8月16日12秒91杰克逊英国1993年8月20日12秒88刘翔中国2006年7月12日12秒87罗伯斯古巴2008年6月12日12秒80梅里特美国2012年9月7日数据：男子110米栏记录问题：估计什么时候突破12秒50？4.1预备知识：最小二乘拟合4.1预备知识：最小二乘拟合假设已知经验公式y=f(c,x)(c为参数,x为自变量),要求根据一批有误差的数据(xi,yi),i=0,1,…,n,确定参数c.这样的问题称为数据拟合。最小二乘法：求c使得平方误差最小化

Q(c)=4.2函数零点MATLAB指令多项式y=polyval(p,x)求得多项式p在x处的值y，x可以是一个或多个点p3=conv(p1,p2)返回多项式p1和p2的乘积[p3,r]=deconv(p1,p2)p3返回多项式p1除以p2的商，r返回余项x=roots(p)求得多项式p的所有复根.p=polyfit(x,y,k)用k次多项式拟合向量数据(x,y),返回多项式的降幂系数MATLAB中一个多项式用系数降幂排列向量来表示。例2.用2次多项式拟合下列数据.

x0.10.20.150-0.20.3

y0.950.840.861.061.500.72M文件eg4_2.m例1.求多项式x3+2x2-5的根

»p=[120-5];x=roots(p),polyval(p,x)Fun=@Mfun定义一个函数句柄，这里Mfun是函数的M文件表达方式Fun=@(var)funstr

定义匿名函数，其中var是变量名,funstr是函数的表达式非线性函数的MATLAB表达

x=fzero(Fun,x0)

返回一元函数Fun的一个零点.x0为标量时,

x

返回函数在x0附近的零点；

x0为向量[a,b]时,

x

返回在[a,b]中的零点(两端函数异号)

4.2函数零点MATLAB指令[x,f,h]=fsolve(Fun,x0)x:

返回多元函数Fun在x0附近的一个零点，其中x,x0均为向量;f:

返回Fun在零点的函数值,应该接近0;h:

返回值如果大于零，说明计算结果可靠，否则计算结果不可靠。例3求函数y=xsin(x2-x-1)在(-2,-0.1)内的零点

>>fun=@(x)x*sin(x^2-x-1)>>fzero(fun,[-2-0.1])>>

fplot(fun,[-2,-0.1]),gridon;>>fzero(fun,[-2,-1.2]),fzero(fun,[-1.2,-0.1])>>fzero(fun,-1.6),fzero(fun,-0.6)例4求方程组在原点附近的解xx(1)yx(2)functionf=eg4_4fun(x)f(1)=4*x(1)-x(2)+exp(x(1))/10-1;f(2)=-x(1)+4*x(2)+x(1)^2/8;>>[x,f,h]=fsolve(@eg4_4fun,[00])使用函数句柄%使用匿名函数fun=@(x)[4*x(1)-x(2)+exp(x(1))/10-1,-x(1)+4*x(2)+x(1)^2/8];[x,f,h]=fsolve(fun,[00])

min(y)

返回向量y的最小值max(y)

返回向量y的最大值[x,f]=fminbnd(fun,a,b)

x返回一元函数fun在[a,b]内的局部极小值点,f返回局部极小值[x,f]=fminsearch(fun,x0)

x返回多元函数fun在初始值x0

附近的局部极小值点,f返回局部极小值.

x,x0均为向量。4.2函数极值MATLAB指令例5.求二元函数f(x,y)=5-x4-y4+4xy在原点附近的极大值。解：maxfmin(-f)x

x(1),y

x(2)>>fun=@(x)x(1)^4+x(2)^4-4*x(1)*x(2)-5;>>[x,g]=fminsearch(fun,[0,0])注：在使用fsolve,fminsearch等指令时，多变量必须合写成一个向量变量，如用x(1),x(2),…。4.2最小二乘拟合MATLAB指令假设已知经验公式y=f(c,x)(c为参数,x为自变量),要求根据一批有误差的数据(xi,yi),i=0,1,…,n,确定参数c.这样的问题称为数据拟合。最小二乘法就是求c使得平方误差最小化

Q(c)=基本格式:c=lsqcurvefit(Fun,c0,x,y)

这里Fun(c,x)为两个输入变量的函数句柄或匿名函数，c0为参数c的预估值，作为迭代初值，x,y为数据向量

完整格式：[c,Q]=lsqcurvefit(Fun,c0,x,y,lb,ub)

lb和ub分别表示c的下界和上界。c返回参数值，Q返回误差平方和。自变量x可以是多变量，这时第三输入参数x应为矩阵。

最小二乘拟合函数类型的选择根据物理或数学机理作图看趋势常用直线、二次函数，指数函数等参数初始值的选择根据参数的实际意义由数据估计网格化搜索常用原点男子110栏问题Matlab程序clear;closeall;data=[13.2413.2113.1613.0012.9312.9212.9112.8812.8712.80];year=[1972197719791979198119891993200620082012];plot(year,data,'o');y=data;x=year-1972;fun=@(c,x)c(1)*exp(-c(2)*x);c=lsqcurvefit(fun,[13,(13.24-12.80)/40],x,y)fun2=@(x)c(1)*exp(-c(2)*x)-12.5;x2=fsolve(fun2,50)%初估50年后，

%计算得73年以后，即2045年xx=0:75;yy=fun(c,xx);figure;plot(year,data,'o',xx+1972,yy);gridon;男子110栏问题的求解模型估计大约2045年突破12秒50

迭代法是从解的初始近似值x0(简称初值)开始，利用某种迭代格式xk+1=g(xk),求得一近似值序列x1,x2,…,xk,xk+1,…逐步逼近于所求的解

(称为不动点)。最常用的迭代法是牛顿迭代法，其迭代格式为

1迭代法4.3计算实验：迭代法Newton法几何意义：切线法切线代替曲线牛顿法程序newton.mfunctionx=newton(fname,dfname,x0,e)ifnargin<4,e=1e-4;end%默认精度x=x0;x0=x+2*e;%使while成立，且进入while后x0得到赋值whileabs(x0-x)>ex0=x;x=x0-fname(x0)/dfname(x0);end例6

求方程x2-3x+ex=2的正根(要求精度

=10-6)解令f(x)=x2-3x+ex-2,f(0)=-1,f(2)>0,所以根在[0,2]内。M文件eg4_6先用图解法找初值，再用牛顿法程序newton.m求解。2．线性化拟合例7.用函数y=aebx

拟合例2的数据方法一：非线性拟合，记a=c(1),b=c(2)>>fun=@(c,x)c(1)*exp(c(2)*x;>>x=…;y=…;[c,Q]=lsqcurvefit(fun,[0,0],x,y)方法二：线性化拟合，两边取对数得z=lny=lna+bx转化为线性拟合。M文件eg4_7不难算出，你向银行总共借了25.2万，30年内共要还51.696万，这个案例中贷款年利率是多少呢？

例8.下面是《新民晚报》2000年3月30日上的一则房产广告：4.4建模实验：购房贷款的利率解设xk为第k个月的欠款数，

a为月还款数,r为月利率。xk+1=(1+r)xk-a那么

xk=(1+r)xk-1-a=(1+r)2xk-2–(1+r)a–a=……=(1+r)kx0

–a[1+(1+r)+……+(1+r)k-1]=(1+r)kx0

–a[(1+r)k-1]/r根据a=0.1436,x0=25.2,x360=0得到25.2(1+r)360

–0.1436[(1+r)360-1]/r=0很难用roots求解！常识上，r应比当时活期存款月利率略高一些。我们用活期存款月利率0.0198/12作为迭代初值，用fzero求解>>clear;fun=@(r)25.2*(1+r)^360-((1+r)^360-…1)/r*0.1436;>>r=fzero(fun,0.0198/12);>>R=12*r得年利率为5.53%.(你知道最新利率吗？)分期付款月还款公式方程(1+r)Nx0–a[(1+r)N-1]/r=0月还款计算x0剩余借款额；N剩余月数；r月利率=年利率/12线性迭代xk+1=axk+b

收敛：当|a|<1,收敛于它的不动点x=b/(1-a)，不收敛：当|a|>1,趋于无穷大。

4.4建模实验：混沌不收敛的非线性迭代可能会趋于无穷大，趋于一个周期解，在一个有限区域内杂乱无章地游荡，这类由确定性运动导致的貌似随机的现象称为混沌现象

4.4建模实验：混沌例（混沌）：昆虫数量的Logistic模型xk+1=axk(1-xk), 0

a

4xk表示第k代昆虫数量(1表示理想资源环境最大可能昆虫数量)。a为资源系数.0

a

4保证了xk在区间(0,1)上封闭。当0

a<1,在[0，1]内有一个不动点0，且由|g’(0)|=a<1,可知它是稳定的,

说明资源匮乏时,昆虫趋于消亡;对于Logistic模型解得有两个不动点0和1-1/a当a>1,不动点0不再稳定;当1<a

3，由|g’(1-1/a)|=|2-a|<1可知不动点1-a-1稳定，说明资源适当时，昆虫稳定于一定数量。a>3，出现两个周期2解，且3<a<1+

周期2轨道稳定。迭代开始发生所谓倍周期分岔，从周期2,周期4,…,周期2n,…

直到a

=3.569945672…。

说明a在[3，a

]取值时，昆虫数量呈现规律性振荡。a>a

迭代序列几乎杂乱无章，即所谓混沌。分岔图程序clear;close;a=0:0.01:4;M=length(a);K=1000;x=zeros(K,M);x(1,:)=rand(1,M);form=1:M,fork=1:K-1x(k+1,m)=a(m)*x(k,m)*(1-x(k,m));end,endfork=1:20,plot(a,x(k,:),'.');title(['k=',int2str(k)]);pause(1);end;plot(a,x(900:K,:),'.');holdoff;分岔图*混沌的特征

(i)初值敏感性:

两个任意近的点出发的两条轨迹迟早会分得很开;蝴蝶效应(Lorenz1963)：一只南美热带雨林中的蝴蝶，偶尔扇动几下翅膀，可能在两周后在美国德克萨斯引起一场龙卷风。

(ii)遍历性:

任意点出发的轨迹总会进入

[0,1]内任意小的开区间。例9（蛛网图）我们用蛛网图来显示混沌的遍历性。

yk=axk(1-xk), xk+1=yk蛛网图正好显示迭代计算x0,y0,x1,y1,……的一系列变化过程。

eg4_9.m蛛网图M函数eg4_9.mfunctionf=eg4_9(a,x0,m,n)x=0:0.01:1;y=a*x.*(1-x);plot(x,x,'r',x,y,'r');holdon;clearx,y;x(1)=x0;y(1)=a*x(1)*(1-x(1));x(2)=y(1);ifm<2,plot([x(1),x(1),x(2)],[0,y(1),y(1)]);endfori=2:ny(i)=a*x(i)*(1-x(i));x(i+1)=y(i);ifi>m,plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]);endendholdoff;习题ex1,ex2,ex6,ex8(1)ex9,ex10(1)(3)ex12，ex13ex13图thR10MATLAB数学实验第五章应用微积分第五章应用微积分5.1预备知识：微积分的基本概念5.2数值微积分MATLAB指令5.3计算实验：数值微积分5.4建模实验：奶油蛋糕、作案时间1.极限和连续数列极限:

>0,

N>0，使当n>N时有

xn-a

<

，则函数极限:如果当x

x0时有f(x)

A，则连续:如果当x

x0时，有f(x)

f(x0)

则称f(x)在x0连续。结论：闭区间上连续函数必有最大值(全局极大值)和最小值(全局极小值)

。5.1预备知识：微积分2.微分与导数函数f(x)在点x=x0的导数为当f’(x0)>0,函数在x0点附近是上升的；当f’(x0)<0,函数在x0点附近是下降的；当f’(x0)=0,x0为驻点（很可能是极值点）函数可导连续几何意义：切线的斜率当n=0得微分中值定理

f(x)-f(x0)=f’(

)(x-x0)

其中

是x0与x之间的某个点。Taylor公式:当f(x)在含有x0某个开区间内具有直到n+1阶的导数，Taylor展开的几何解释任意函数的局部多项式逼近3.多元函数微分学

设f(x,y)在点(x0,y0)附近有定义，当(x,y)以任何方式趋向于(x0,y0)时，f(x,y)趋向于一个确定的常数A，则二重极限若A=f(x0,y0),称f(x,y)在(x0,y0)点连续f(x,y)在点(x0,y0)的偏导数多元函数Taylor公式

函数极值一元可微f(x)在内点x0取得局部极值必要条件是f’(x0)=0.局部极大充分条件是f’(x0)=0,f’’(x0)<0局部极小充分条件是f’(x0)=0,f’’(x0)>0.函数极值二元可微f(x,y)在内点(x0,y0)取得局部极大或极小的必要条件是梯度

f(x0,y0)=0,局部极大(或局部极小)充分条件是

f(x0,y0)=0且下列Hesse矩阵负定（或正定）。梯度

f=(),Hesse矩阵多元函数极值也有类似的结果。4.积分

函数f(x)在区间[a,b]上的定积分其中a=x0<x1<…<xn=b,

xi=xi-xi-1,

i

(xi-1,xi),i=1,2,…,n几何意义：曲边梯形的面积牛顿-莱布尼兹公式：若在[a,b]上,F’(x)=f(x),则二重积分定义曲线&曲面平面曲线(x(t),y(t)),a<t<b的长度为空间曲线(x(t),y(t),z(t)),a<t<b的长度为曲面z=z(x,y),(x,y)

G的面积为1.数值差分

n维向量x=(x1,x2,

,xn)的差分定义为n-1维向量

x=(x2-x1,x3-x2,

,xn-xn-1)。diff(x)如果x是向量，返回向量x的差分如果x是矩阵，则按各列作差分。diff(x,k)k阶差分，即差分k次5.2数值微积分MATLAB指令2.数值导数和梯度q=polyder(p)

求得由向量p表示的多项式导函数的向量表示q.Fx=gradient(F,x)

返回向量F表示的一元函数沿x方向的导函数F’(x).

其中x是与F同维数的向量.[Fx,Fy]=gradient(F,x,y)

返回矩阵F表示的二元函数的数值梯度(F’x,F’y),当F为m×n

矩阵时,x,y分别为n维和m维的向量.例子»clear;x=[11.11.21.3];y=x.^3;»3*x.^2»dy=diff(y)./diff(x)»dy=gradient(y,x)注意(1)gradient内点用中心差商精度高，两端用向前或向后差商精度低；(2)计算精度随着步长变小而提高。3.梯形积分法

z=trapz(x,y)

返回积分的近似值，其中x表示积分区间的离散化向量;y是与x同维数的向量,表示被积函数。例1解»clear;x=-1:0.1:1;y=exp(-x.^2);

»trapz(x,y)注意：积分时被积函数要用数组运算!!4.高精度数值积分

z=quadl(Fun,a,b).z=integral(Fun,a,b)不是数字1，是字母！>>z=quadl(@(x)exp(-x.^2),-1,1)>>z=integral(@(x)exp(-x.^2),-1,1)注意：积分中fun里用数组运算!!>>fun=@(x)exp(-x.^2).*log(x).^2;>>q=integral(fun,0,Inf)%integral能求反常积分

5.

重积分z=dblquad(Fun,a,b,c,d)

求得二元函数Fun(x,y)在矩形区域的重积分.z=triplequad(Fun,a,b,c,d,e,f)求得三元函数Fun(x,y,z)在长方体区域上的三重积分。z=quad2d(Fun,a,b,cx,dx)求得二元函数Fun(x,y)的一般区域重积分。a,b为变量x的下、上限；cx,dx为变量y的下、上限函数(自变量为x)z=integral2(Fun,a,b,cx,dx)类似quad2dz=integral3(Fun,a,b,cx,dx,exy,fxy)求得三元函数Fun(x,y,z)一般区域的三重积分例5.2计算重积分fun=@(z,x,y)log(3+x+y+z)+z.^2.*cos(x);xmin=@(z)-sqrt(1-z.^2);xmax=@(z)sqrt(1-z.^2);ymin=@(z,x)-sqrt(1-x.^2-z.^2);ymax=@(z,x)sqrt(1-x.^2-z.^2);q=integral3(fun,0,1,xmin,xmax,ymin,ymax)fun=@(x,y)1./(sqrt(x+y).*(1+x+y).^2);q=integral2(fun,0,1,0,@(x)1-x)1.数值微分若f(x)在x=a可导,设h>0且足够小称为向前差商向后差商中心差商5.3计算实验：数值微积分差商公式比较中心差商精度比较高自动步长的中心差商程序：deriv.mfunctiond=deriv(fname,a,h0,e)h=h0;d=(fname(a+h)-fname

(a-h))/2/h;d0=d+2*e;%为了保证abs(d-d0)>e成立whileabs(d-d0)>ed0=d;h0=h;h=h0/2d=(fname(a+h)-fname(a-h))/2/h;end>>deriv(@(x)sin(x),pi/6,0.1,1e-4)例5.4（奶油蛋糕）某数学家的学生要送一个特大的蛋糕来庆贺他90岁生日。为了纪念他提出的口腔医学的悬链线模型，学生们要求蛋糕店老板将蛋糕边缘半径作成下列悬链线函数

r=2-(exp(2h)+exp(-2h))/5,0<h<1(单位：米)。问如何计算重量？5.4建模实验：解设高为H，半径r,比重为k若蛋糕是单层圆盘的，则蛋糕的重量为：

W=k

Hr2rHr1r2若蛋糕是双层的，每层高H/2，下层半径r1,上层半径r2，则

W=k

H(r12+r22)/2如果蛋糕是n层的，每层高H/n，半径分别r1,…,rn,则fun=@(h)(2-(exp(2*h)+exp(-2*h))/5);n=5;h=0.1:0.2:0.9;r=fun(h);W=pi*1/n*sum(r.^2)计算结果5.4327当n

，fun2=@(h)(2-(exp(2*h)+exp(-2*h))/5).^2;pi*integral(fun2,0,1)计算结果5.4171即旋转体体积公式解：t=0时z=R.考虑在时间dt内水面变化dz，漏水的体积为

uAdt=-

x2dz其中u为出水速度，A为小孔面积，x为高度z水面的半径5m0.5m0例5.5（作案时间）锅炉房中发生一起命案，球形锅炉底部有个弹孔在漏水。锅炉半径R=5m，平时充满了水，测得底部小孔半径b=0.1m，此时水面已将下降至离底部0.5m。推测枪击是多长时间以前发生的?分析：求时间t与水面高度z的对应关系。z表示从球心测量的水面高度A=

b2,x2=R2-z2注意：水从孔漏出的速度u不是常数！zxt=0时z=R.

在顶部水降到0.5m时，z=0.5-R，从而

t=0+5m0.5m0u由下列能量方程决定mg(z+R)=mu2/2

m是质量，g为重力加速度。%M脚本eg5_6.mclear;R=5;b=0.1;g=9.81;z1=0.5-R;z2=R;n=100;h=(z2-z1)/n;z=z1:h:z2;f=(R^2-z.^2)./(b^2*sqrt(2*g*(z+R)));I=trapz(z,f)/60/60结果为0.5144小时。Page101习题ex4ex5(2)(4)(6)(7)(8)ex6ex10ex12ex13MATLAB数学实验

第六章常微分方程第六章常微分方程6.1预备知识：常微分方程6.2解常微分方程的MATLAB指令6.3计算实验：Euler法和刚性方程组6.4建模实验:导弹系统的改进1.微分方程的概念常微分方程:f(t,y,y’,y’’,…,y(n))=0常微分方程的解:函数y(t).6.1预备知识：常微分方程微分方程组:联系一些未知函数x(t),y(t),z(t),…的一组微分方程.偏微分方程:含有偏导数的微分方程.其解为多元函数u(t,x,y,z)。6.1预备知识：常微分方程2.常微分方程解析求解：某些特殊的方程(1)初等积分法例如y’=ty,y(0)=1.分离变量dy/y=tdt,lny=t2/2+c,由初始值得c=0,所以y(t)=exp(t2/2)（2）线性常系数微分方程的特征根法线性方程:y(n)+a1(t)y(n-1)+…+an-1(t)y’+an(t)y=b(t)常系数方程:

若ai(t)(i=1,…,n)与t无关。齐次方程:

若b(t)=0。

y(n)+a1y(n-1)+…+an