清华大学数学实验报告6
数学建模基础实验报告(3篇)
第1篇一、实验目的本次实验旨在让学生掌握数学建模的基本步骤，学会运用数学知识分析和解决实际问题。
通过本次实验，培养学生主动探索、努力进取的学风，增强学生的应用意识和创新能力，为今后从事科研工作打下初步的基础。
二、实验内容本次实验选取了一道实际问题进行建模与分析，具体如下：题目：某公司想用全行业的销售额作为自变量来预测公司的销售量。
表中给出了1977—1981年公司的销售额和行业销售额的分季度数据（单位：百万元）。
1. 数据准备：将数据整理成表格形式，并输入到计算机中。
2. 数据分析：观察数据分布情况，初步判断是否适合使用线性回归模型进行拟合。
3. 模型建立：利用统计软件（如MATLAB、SPSS等）进行线性回归分析，建立公司销售额对全行业的回归模型。
4. 模型检验：对模型进行检验，包括残差分析、DW检验等，以判断模型的拟合效果。
5. 结果分析：分析模型的拟合效果，并对公司销售量的预测进行评估。
三、实验步骤1. 数据准备将数据整理成表格形式，包括年份、季度、公司销售额和行业销售额。
将数据输入到计算机中，为后续分析做准备。
2. 数据分析观察数据分布情况，绘制散点图，初步判断是否适合使用线性回归模型进行拟合。
3. 模型建立利用统计软件进行线性回归分析，建立公司销售额对全行业的回归模型。
具体步骤如下：（1）选择合适的统计软件，如MATLAB。
（2）输入数据，进行数据预处理。
（3）编写线性回归分析程序，计算回归系数。
（4）输出回归系数、截距等参数。
4. 模型检验对模型进行检验，包括残差分析、DW检验等。
（1）残差分析：计算残差，绘制残差图，观察残差的分布情况。
（2）DW检验：计算DW值，判断随机误差项是否存在自相关性。
5. 结果分析分析模型的拟合效果，并对公司销售量的预测进行评估。
四、实验结果与分析1. 数据分析通过绘制散点图，观察数据分布情况，初步判断数据适合使用线性回归模型进行拟合。
2. 模型建立利用MATLAB进行线性回归分析，得到回归模型如下：公司销售额 = 0.9656 行业销售额 + 0.01143. 模型检验（1）残差分析：绘制残差图，观察残差的分布情况，发现残差基本呈随机分布，说明模型拟合效果较好。
高等数学数学实验报告(两篇)
引言概述:高等数学数学实验报告（二）旨在对高等数学的相关实验进行探究与研究。
本次实验报告共分为五个大点，每个大点讨论了不同的实验内容。
在每个大点下，我们进一步细分了五到九个小点，对实验过程、数据收集、数据分析等进行了详细描述。
通过本次实验，我们可以更好地理解高等数学的概念和应用。
正文内容:一、微分方程实验1.利用欧拉法求解微分方程a.介绍欧拉法的原理和步骤b.详细阐述欧拉法在实际问题中的应用c.给出具体的实例，展示欧拉法的计算步骤2.应用微分方程建立模型求解实际问题a.介绍微分方程模型的建立方法b.给出一个具体的实际问题，使用微分方程建立模型c.详细阐述模型求解步骤和结果分析3.使用MATLAB求解微分方程a.MATLAB求解微分方程的基本语法和函数b.给出一个具体的微分方程问题，在MATLAB中进行求解c.分析结果的准确性和稳定性二、级数实验1.了解级数的概念和性质a.简要介绍级数的定义和基本概念b.阐述级数收敛和发散的判别法c.讨论级数的性质和重要定理2.使用级数展开函数a.介绍级数展开函数的原理和步骤b.给出一个函数，使用级数展开进行近似计算c.分析级数近似计算的精确度和效果3.级数的收敛性与运算a.讨论级数收敛性的判别法b.介绍级数的运算性质和求和法则c.给出具体的例题，进行级数的运算和求和三、多元函数极值与最值实验1.多元函数的极值点求解a.介绍多元函数的极值点的定义和求解方法b.给出一个多元函数的实例，详细阐述求解过程c.分析极值点对应的函数值和意义2.多元函数的条件极值与最值a.讨论多元函数的条件极值的判定法b.给出一个具体的多元函数，求解其条件极值和最值c.分析条件极值和最值对应的函数值和意义3.利用MATLAB进行多元函数极值与最值的计算a.MATLAB求解多元函数极值与最值的基本语法和函数b.给出一个多元函数的具体问题，在MATLAB中进行求解c.分析结果的准确性和可行性四、曲线积分与曲面积分实验1.曲线积分的计算方法与应用a.介绍曲线积分的定义和计算方法b.给出一个具体的曲线积分问题，详细阐述计算过程c.分析曲线积分结果的几何意义2.曲线积分的应用举例a.讨论曲线积分在实际问题中的应用b.给出一个实际问题，使用曲线积分进行求解c.分析曲线积分结果的实际意义和应用价值3.曲面积分的计算方法与应用a.介绍曲面积分的定义和计算方法b.给出一个具体的曲面积分问题，详细阐述计算过程c.分析曲面积分结果的几何意义五、空间解析几何实验1.空间曲线的参数方程表示与性质a.介绍空间曲线的参数方程表示和性质b.给出一个具体的空间曲线，转化为参数方程表示c.分析参数方程对应的几何意义和性质2.平面与空间直线的位置关系a.讨论平面与空间直线的位置关系的判定方法b.给出一个具体的平面与空间直线的问题，判定其位置关系c.分析位置关系对应的几何意义和应用实例3.空间直线与平面的夹角和距离计算a.介绍空间直线与平面的夹角和距离的计算方法b.给出一个具体的空间直线和平面，计算其夹角和距离c.分析夹角和距离计算结果的几何意义总结:通过本次高等数学数学实验报告（二），我们深入了解了微分方程、级数、多元函数极值与最值、曲线积分、曲面积分以及空间解析几何的相关概念和应用。
回归分析-数学实验-清华大学
回归分析2012011849 分2 李上【实验目的】1.了解回归分析的基本原理，掌握MATLAB实现的方法；2.练习用回归分析解决实际问题。
【实验内容】1.题目5问题描述社会学家认为犯罪与收入低、失业及人口规模有关，对20个城市的犯罪率y（每10万人中犯罪的人数）与年收入低于5000美元家庭的百分比x1、失业率x2和人口总数x3（千人）进行了调查（表格略）（1）若x1~x3中至多只许选择2个变量，最好的模型是什么？（2）包含3个自变量的模型比上面的模型好吗？确定最终模型。
（3）对最终模型观察残差，有无异常点，若有，剔除后如何。
问题分析先做y和x i的散点图，来大致判断自变量和因变量的关系。
而后选择x i中的某些与y进行线性回归并与包含三个自变量的模型进行比较，并剔除某些数据点，得出更好的拟合结果。
代码和结果第一问y=[11.2 14.5 12.7 28.9 13.4 26.9 20.9 14.9 40.7 15.7 35.7 25.8 5.3 36.2 8.7 21.7 24.8 18.1 9.6 25.7]x1=[16.5 18.1 16.5 24.9 20.5 23.1 20.2 17.9 26.3 19.1 21.3 22.4 16.5 24.7 17.2 20.2 19.2 18.6 14.3 16.9] x2=[6.2 6 5.9 8.3 6.4 7.4 6.4 6.7 9.3 5.8 7.6 8.6 5.3 8.6 4.9 8.4 7.3 6.5 6.4 6.7]x3=[587 7895 643 854 643 762 1964 716 635 2793 1531921 692 741 713 595 1248 625 49 3353]plot(y,x1,'r+');pauseplot(y,x2,'r+');pauseplot(y,x3,'r+');图像如下：y-x1图像y-x2图像y-x3图像因此，y和x1、x2的关系大致为线性关系，所以，选择x1和x2和y做二元线性回归。
数学建模实验报告实验六
《国土面积计算的研究方案》实验报告实验项目：插值与拟合问题实验地点：实验室名称：学院：管理科学与工程学院年级专业班：工程管理111学生姓名：刘继壮学号： ********完成时间：2013-04-20 教师评语：开课时间：至学年第学期成绩教师签名批阅日期国土面积计算的研究方案1 问题重述已知欧洲一个国家的地图，为了算出它的国土面积和边界长度，首先对地图作如下测量：以由西向东方向为x 轴，由南向北方向为y 轴，选择方便的原点，并将从最西边界点到最东边界点在x 轴上的区间适当地分为若干段，在每个分点的y 方向测出南边界点和北边界点的y 坐标1y 和2y ，这样就得到了表8的数据（单位：mm ）。
根据地图的比例我们知道1mm 相当于4km 。
请你研究如下三个问题：（1）直接进行数值积分求该国边界的近似长度和国土的近似面积。
（2）先对边界进行三次样条插值，再计算该国边界的近似长度和国土的近 似面积。
（3）对（1）和（2）的计算结果进行比较。
2 模型的建立与求解解： 假设测量的地图和数据准确，由最西边界点与最东边界点分为上下两条连续的边界曲线，边界内的所有土地均为该国国土。
假设从最西边界点到最东边界点，变量[]b a x ,∈，划分[]b a x ,∈为n 小段[]i i x x ,1-，并由此将国土分为n 小块，设每一小块均为X型区域。
即做垂直于x 轴的直线穿过该区域，直线与边界曲线最多只有两个交点。
利用数值积分法将上边界点与下边界点分别利用插值函数求出两条曲线，则曲线所围成面积即为国土面积（地图上的国土面积），然后根据比例缩放关系求出国土面积的近似解。
设上边界函数为()x f 2，下边界函数为()x f 1，由定积分定义可知曲线所围区域面积为()()ini i i n b a x f f dx x f ∆-=∑⎰=+∞→112)]([lim εε式中，],[1i i i x x -∈ε利用软件求得的图形为：020*********12014016020406080100120140东西距离(单位:mm)南北距离(单位:m m )国土面积计算——三次插值(比例尺为1:4000000)某国家国土区域插值边界参考文献[1]《运筹学》教材编写组，运筹学（修订版），北京：清华大学出版社，1990。
实验报告——线性规划建模与求解
exitflag =1
实验过程记录（含：基本步骤、主要程序清单及异常情况记录等）（接上页）：
实验书中的实际问题求解：
解：设a 为0-1变量，表示第i根8M线材
设b 为0-1变量，表示第i根12M线材
X 表示第i根8M线材截得的第j种长度的线材数目
Y 表示第i根12M线材截得的第j种长度的线材数目
5.完成实验中的实际问题求解。
实验过程记录（含：基本步骤、主要程序清单及异常情况记录等）：
习题求解
1.2将下列线性规划转化为标准型，并用程序求解。
解：转化为标准型如下：
用matlab求解命令如下：
f=[-3,4,-2,5,0,0];
aeq=[4,-1,2,-4,0,0;1,1,2,-1,1,0;-2,3,-1,2,0,-1];
b=[-60,-70,-60,-50,-20,-30]’;
lb=zeros(6,1);
[x,fval,exitflag,output,lambda]=linprog(f,a,b,[],[],lb);
解得结果为：
x =[41.9176,28.0824,35.0494,14.9506,9.8606,20.1394]
Z为浪费的线材总长度
又由于150*(8+12)远大于所需线材总长度，故知所用两种线材每种不超过150根
解不出
实验结果报告与实验总结：
对于实验指导书中matlab使用的例题和方法已经基本掌握，《运筹学》书中例题与方法处于基本了解的程度，不能灵活运用，但书后习题全都能独立完成，已经有一定解题能力。且实验书中的实际运用题的简易版问题的解题方法也已经掌握，但此实验题仍很吃力。
fval = 3.6000
填料塔实验报告 清华大学
填料吸收塔实验报告朱奇化22 2012011910同组成员时羽剑实验日期2014-05-10 实验目的1.了解填料吸收塔的构造并实际操作；2.了解填料塔的流体力学性能。
实验原理压强降是塔设计中的重要参数，气体通过填料层压强降的大小决定了塔的动力消耗。
压强降与气液流量有关，不同喷淋量下填料层的压强降ΔP与空塔气速u的关系如下图所示：图6-1 填料层的ΔP～u关系当无液体喷淋即喷淋量L=0时，干填料的ΔP～u的关系是直线，如图中的直线0。
当有一定的喷淋量时，ΔP～u的关系变成折线，并存在两个转折点，下转折点称为“载点”，上转折点称为“泛点”。
这两个转折点将ΔP～u关系分为三个区段：恒持液量区、载液区与液泛区。
实验步骤1.测量干填料层(△P／Z)─u关系曲线:先全开调节阀2,后启动鼓风机,用阀2调节进塔的空气流量,按空气流量从小到大的顺序读取填料层压降△P,转子流量计读数和流量计处空气温度,测量12～15组数据•然后在双对数坐标纸上以空塔气速u为横坐标,以单位高度的压降△P／Z为纵坐标,标绘干填料层(△P／Z)─u关系曲线。
2.测量某喷淋量下填料层(△P／Z)─u关系曲线:用水喷淋量为600L／h时,用上面相同方法读取填料层压降△P,•转子流量计读数和流量计处空气温度并注意观察塔内的操作现象, •一旦看到液泛现象时记下对应的空气转子流量计读数。
在对数坐标纸上标出液体喷淋量为600L／h下(△P／z)─u•关系曲线,确定液泛气速并与观察的液泛气速相比较。
3.测量某喷淋量下填料层(△P／Z)─u关系曲线:用水喷淋量为1000L／h时,用上面相同方法读取填料层压降△P,•转子流量计读数和流量计处空气温度并注意观察塔内的操作现象, •一旦看到液泛现象时记下对应的空气转子流量计读数。
在对数坐标纸上标出液体喷淋量为1000L／h下(△P／z)─u•关系曲线,确定液泛气速并与观察的液泛气速相比较。
4.实验完毕后，关闭空压机、真空泵、进水阀门、等仪器设备的电源，并将所有仪器复原。
大学数学实验课后习题答案(清华大学出版)
大学数学实验课后习题答案(清华大学出版)实验名称：MA TLAB 程序设计(1)作马鞍面：22,66,8823x y z x y =--≤≤-≤≤程序: x=-6:0.5:6;y=-8:0.5:8[X,Y]=meshgrid(x,y);Z=X.^2./2-Y .^2./3;mesh(X,Y ,Z)(2)P441第5题程序1:n=18;I(1)=1-exp(-1);%I(1)对应I0for k=1:n-1I(k+1)=1-(k+1)*I(k);end I程序2:n=18;I1=(1/(n+1))*exp(-1);I2=1/(n+1);I(18)=(I1+I2)/2;for k=n:-1:2I(k-1)=(1-I(k))/n;endI（3）自定义函数：lnsin cos ln tan y x x x =-，并求()?3y π =程序：function y=fun(x);y=log(sin(x))-cos(x)*log(tan(x));>>fun(pi/3)（4）P441第10题的（1）、（2）小题。
要求建立函数M 文件求解。
并求：201!n T n ==∑程序1：求!n 自定义函数function y=fun(n)A=1;for k=1:nA=A*k;endA程序2：求：201!n T n ==∑s=0; for n=1:20A=1;for k=1:nA=A*k;ends=s+A;endsC程序3：求nmfunction y=funa(n,m) A=1;%求for k=1:nA=A*k;endB=1;for k=1:mB=B*k;endC=1;for k=1:n-mC=C*k;endD=A/(B*C) %求组合数一元函数的图形练习解答： 1．用ezplot画出的图象.程序：ezplot('asin(x)') 2．用ezplot画出用在(0,)之间的图象.程序：ezplot('sec(x)',[0 pi])3．在同一坐标系中画出,,,,的图象.并用gtext加以标记ezplot('sqrt(x)')hold onezplot('x^2')hold onezplot('x^(1/3)')hold onezplot('x^3')hold onezplot('x')axis([-2 3 -2 2])gtext('sqrt(x)')gtext('x^2')gtext('x^(1/3)')gtext('x^3')gtext('x')4．画出及其反函数的图象. x=-2:0.01:20;y=1+log(x+2+eps);plot(x,y)holdplot(y,x,'r')axis([-4 4 -4 4])8题：x=100;y=50;n=50; r1=0.2;r2=0.3;a1=0.001;a2=0.002;for k=1:nx(k+1)=(1+r1-a1*y(k))*x(k);y(k+1)=(1-r2+a2*x(k))*y(k);endk=0:n;round([k',x',y'])plot(k,x,k,y),grid,2题：function z=exf14(x0,y0,n,r,N,d,a,b); x=x0;y=y0;for k=1:nx(k+1)=x(k)+r*(1-x(k)/N)*x(k)-a*y(k)*x(k)/N; y(k+1)=(1-d+b*x(k)/N)*y(k);endz=[x',y'];z=exf14(1000,100,100,0.8,3000,0.9,1.6,1.5); k=0:100;plot(k,z(:,1),k,z(:,2)),grid。
清华大学偏振光学实验完整实验报告
偏振光学实验完整实验报告工物53 李哲 2015011783 16号1.实验目的：（1）理解偏振光的基本概念，在概念以及原理上了解线偏振光，圆偏振光以及椭圆偏振光，并了解偏振光的起偏与检偏方法。
以及线偏振光具有的一些性质。
（2）学习偏振片与玻片的工作原理。
2.实验原理：（1）光波偏振态的描述：· 单色偏振光可以分解成两个偏振方向垂直的线偏振光的叠加：t a E X ωcos 1=与()δω+=t a E Y cos 1（其中δ是两个偏振方向分量的相位延迟，21,a a 为两个光的振幅），由其中的δ,,21a a 就可以确定这个线偏振光的性质。
πδ=或0=δ就为线偏振光，2,21πδ==a a 为圆偏振光（就是光矢量的顶点绕其中点做圆周运动，依然是偏振光），而一般情况下是椭圆偏振光。
· 上述式子通常描述的是椭圆偏振光，而本实验通过测量椭圆的长轴方位角ψ以及椭圆的短半轴与长半轴的比值对于椭圆偏振光进行描述。
其计算式是：()δβcos 2tan arctan 21⋅=ψ()12sin sin 112222-⋅-+=βδa b而对于实验中的椭圆偏振光而言，其光强在短轴对应的方向最小，在长轴的对应方向最大，所以可以通过使这个椭圆偏振光通过一个偏振片，并调整偏振片的透射轴方位，测量其最大最小值，就可以知道其长轴短轴的比值。
又由于光强与振幅的平方成正比，所以测得的光强的比值是长轴短轴之比的平方。
（2）偏振片：· 理想偏振片：只有电矢量振动方向与透射轴平行方向的光波分量才能通过偏振片。
· 实验中的偏振片不是理想化的，并不能达到上述的效果，当入射光波的振动方向与透射轴平行时，其透射率不能达到1，当垂直于透射轴时，其透射率不是0。
所以对于偏振片有主透射率以及消光比两个量进行描述。
· 主透射率21T T ，指沿透射轴或消光轴方向振动光的光强透射率。
两者的比值是消光比e 。
清华大学 数学实验06讲
Experiments in Mathematics
实验6 非线性方程求解 清华大学数学科学系
什么叫方程（组）?
方程：包含未知数（或/和未知函数）的等式 方程组：包含未知数（或/和未知函数）的一组等式
不包含未知函数的方程组的一般形式： F(x)=0 x= (x1,x2,…,xn)T, F(x)=(f1(x),f2(x),…,fm(x))T (一般m=n)
f ( xk −1 )
y=f(x)
收敛速度比牛顿切线法稍慢
O xk-1 xk
xk+1
Q
x x*为单根时收敛阶数是 1.618
解非线性方程组的牛顿法
解方程 f(x)=0 的牛顿切线法 推广到解方程组
f (xk+1) = f (xk ) + f ′(xk )(xk +1 − xk ) xk +1 = xk − f (xk ) f ′(xk )
f (3) = −2, f (4) = 6 存在根 x ∈ (3,4)
x
= ϕ (x) 1
= 14− x2 ,
迭代公式：
x k
+1
= 14 − x2 k
x =ϕ (x) =14/(x +1)，迭代公式：x =14/(x +1)
2
k +1
k
x = ϕ (x) = x − (x 2 + x − 14) /(2x + 1) , 3
i =1
非线性方程的基本解法
解方程 f(x)=0第一步——确定根的大致范围
• 图形法：作f(x)图形,观察f(x)与x轴的交点
• 根的隔离:二分法 图形法
x6 − 2x4 − 6x3 −13x2 + 8x +12 = 0
微分方程数值解实验报告
微分方程数值解实验报告微分方程数值解实验报告一、引言微分方程是数学中一类重要的方程，广泛应用于各个科学领域。
在实际问题中，往往难以得到微分方程的解析解，因此需要借助数值方法来求解。
本实验旨在通过数值解法，探索微分方程的数值解及其应用。
二、数值解法介绍常用的微分方程数值解法有欧拉法、改进欧拉法、四阶龙格-库塔法等。
在本实验中，我们将采用改进欧拉法进行数值解的求取。
改进欧拉法是一种一阶的显式迭代法，其基本思想是将微分方程的导数用差商来近似表示，并通过迭代逼近真实解。
具体迭代公式如下：\[y_{n+1} = y_n + h \cdot f(x_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot f(x_n, y_n))\]其中，\(y_n\)表示第n步的近似解，\(h\)表示步长，\(f(x_n, y_n)\)表示微分方程的导数。
三、实验步骤1. 确定微分方程及初始条件在本实验中，我们选择经典的一阶常微分方程：\[y' = -2xy\]并给定初始条件\(y(0) = 1\)。
2. 设置步长和迭代次数为了得到较为准确的数值解，我们需要合理选择步长和迭代次数。
在本实验中，我们将步长设置为0.1，迭代次数为10。
3. 迭代计算数值解根据改进欧拉法的迭代公式，我们可以通过编写计算程序来求解微分方程的数值解。
具体计算过程如下：- 初始化：设定初始条件\(y_0 = 1\)，并给定步长\(h = 0.1\)。
- 迭代计算：使用改进欧拉法的迭代公式，通过循环计算得到\(y_1, y_2, ...,y_{10}\)。
- 输出结果：将计算得到的数值解输出，并进行可视化展示。
四、实验结果与分析通过以上步骤，我们得到了微分方程的数值解\(y_1, y_2, ..., y_{10}\)。
将这些数值解进行可视化展示，可以更直观地观察到解的变化趋势。
根据实验结果，我们可以发现随着迭代次数的增加，数值解逐渐逼近了真实解。
混沌系统实验报告
一、实验目的1. 了解混沌现象的基本概念和特性。
2. 掌握混沌系统实验的基本方法和步骤。
3. 通过实验观察混沌现象，验证混沌系统的基本特性。
4. 理解混沌现象在实际应用中的意义。
二、实验原理混沌现象是自然界和人类社会普遍存在的一种复杂现象，具有以下基本特性：1. 敏感性：初始条件的微小差异会导致系统行为的巨大差异。
2. 无序性：混沌系统表现出复杂、不规则的行为，难以预测。
3. 非线性：混沌系统内部存在非线性相互作用，导致系统行为复杂。
4. 吸引子：混沌系统最终会收敛到一个或多个吸引子上，形成稳定的动态行为。
本实验主要研究一个典型的混沌系统——洛伦茨系统，其数学模型如下：\[\begin{cases}\frac{dx}{dt} = \sigma(y - x) \\\frac{dy}{dt} = x(\rho - z) - y \\\frac{dz}{dt} = xy - \beta z\end{cases}\]其中，\(x\)、\(y\)、\(z\) 分别代表洛伦茨系统的三个状态变量，\(\sigma\)、\(\rho\)、\(\beta\) 为系统参数。
三、实验仪器与设备1. 混沌系统实验仪2. 数字示波器3. 计算机及数据采集软件四、实验步骤1. 打开混沌系统实验仪，连接好实验仪器。
2. 设置洛伦茨系统的参数，包括 \(\sigma\)、\(\rho\)、\(\beta\)。
3. 通过实验仪观察洛伦茨系统的动态行为，并记录实验数据。
4. 使用数字示波器观察洛伦茨系统的相图和时序图。
5. 使用数据采集软件记录洛伦茨系统的状态变量随时间的变化曲线。
6. 分析实验数据，验证混沌系统的基本特性。
五、实验结果与分析1. 当 \(\sigma = 10\)、\(\rho = 28\)、\(\beta = 8/3\) 时，洛伦茨系统呈现出典型的混沌现象。
从时序图可以看出，系统状态变量 \(x\)、\(y\)、\(z\) 随时间的变化呈现出无规则、复杂的振荡行为。
清华大学MATLAB数学实验作业参考
数学实验内容及目的(作业参考)实验1 一元函数的图形实验目的１．学习 matlab一元函数绘图命令.２．进一步理解函数概念.实验内容１．学习matlab命令.matlab绘图命令比较多，我们选编一些常用命令，并简单说明其作用，这些命令的调用格式，可参阅例题及使用帮助help查找.表1.1 二维绘图函数表1.2 基本线型和颜色表1.3 二维绘图工具表1.4 axis 命令linspace 创建数组命令，调用格式为:x=linspace(x1，x2，n)，创建了x 1到x2之间有n 个数据的数组.funtool 函数工具，在matlab 指令窗键入funtool 可打开“函数计算器”图形用户界面.２．绘制函数图形举例. 例1.1.画出x y sin =的图形解:首先建立点的坐标，然后用plot 命令将这些点绘出并用直线连接起来，采用中学五点作图法，选取五点)0,0(、)1,2(π、)0,(π、)1,23(-π、)0,2(π.输入命令:x=[0，pi/2，pi ，3*pi/2，2*pi]；y=sin(x)；plot(x ，y) 这里分号表示该命令执行结果不显示.可以想象，随点数增加，图形越来越接近x y sin =的图象.例如，在0到π2之间取30个数据点，绘出的图形与x y sin =的图象已经非常接近了． x=linspace(0，2*pi ，30)；y=sin(x)；plot(x ，y) 也可以如下建立该图形.x=0:0.1:2*pi ；y=sin(x)；plot(x ，y) 还可以给图形加标记、格栅线. x=0:0.1:2*pi ； y=sin(x)；plot(x ，y ，’r—‘) title(‘正弦曲线’) xlabel(‘自变量 x’) ylabel(‘函数y=sinx’) text(5.5，0，’y=sinx’) grid 上述命令第三行选择了红色虚线，第四行给图加标题“正弦曲线”，第五行给x 轴加标题“自变量x ”，第六行给y 轴加标题“函数x y sin =”，第七行在点)0,5.5(处放置文本“x y sin =”，第八行给图形加格栅线.例1.2.画出xy 2=和xy )2/1(=的图象.解:输入命令:x=-4:0.1:4；y1=2.^x ；y2=(1/2).^x ；plot(x ，y1，x ，y2)； axis([-4，4，0，8])matlab 允许在一个图形中画多条曲线.plot(x1，y1，x2，y2，x3，y3)指令绘制)2(2),1(1x f y x f y ==等多条曲线.matlab 自动给这些曲线以不同颜色.例1.3.画出arctgx y =的图象. 解: 输入命令:x=-20:0.1:20；y=atan(x)；plot(x ，y ，[-20，20]，[pi/2，pi/2]，[-20，20]，[-pi/2，-pi/2]) grid从图上看，arctgx y =是有界函数，2π±=y 是其水平渐近线. 例1.4.在同一坐标系中画出tgx y x y x y ===,,sin 的图象.解:输入命令:x=-pi/2:0.1:pi/2；y1=sin(x)；y2=tan(x)； plot(x ，x ，x ，y1，x ，y2) axis equalaxis([-pi/2，pi/2，-3，3]) grid从图上看，当0>x 时，tgx x x <<sin ，当0<x 时，tgx x x >>sin ，x y =是x y sin =和tgx y =在原点的切线，因此，当1<<x 时，x tgx x x ≈≈,sin .例1.5.画出110-=xy 及)1lg(+=x y 的图形.解:输入命令:x1=-1:0.1:2；y1=10.^x1-1；x2=-0.99:0.1:2；y2=log10(x2+1)； plot(x1，y1，x2，y2)从图上看，这两条曲线与我们所知的图象相差很远，这是因为坐标轴长度单位不一样的缘故.110-=xy 与)1lg(+=x y 互为反函数，图象关于x y =对称，为更清楚看出这一点，我们再画出x y =的图象. hold onx=-1:0.01:2；y=x ；plot(x ，y ，’r’) axis([-1，2，-1，2]) axis square ；hold offplot 语句清除当前图形并绘出新图形，hold on 语句保持当前图形.例1.6.画出心形线)cos 1(3a r +=的图象. 解:输入命令:x=-2*pi:0.1:2*pi ；r=3*(1+cos(x))；polar(x ，r)例1.7.画出星形线t y t x 33sin 3,cos 3==的图象. 解:这是参数方程，可化为极坐标方程.233232)sin (cos 3a a r +=输入命令:x=0:0.01:2*pi ；r=3./(((cos(x)).^2).^(1/3)+((sin(x)).^2).^(1/3)).^(3/2)； polar(x ，r)注意，如果建立r 与t 的关系，此时t 只是参数，不是极坐标系下的极角.实验2 函数极限实验目的１．理解极限概念.２．掌握用matlab 软件求函数极限的方法.实验内容１．学习matlab 命令.matlab 求极限命令可列表如下:表2.1matlab 代数方程求解命令solve 调用格式. solve (函数)(x f ) 给出0)(=x f 的根.２．理解极限概念.数列}{n x 收敛或有极限是指当n 无限增大时，n x 与某常数无限接近或n x 趋向于某一定值，就图形而言，也就是其点列以某一平行与y 轴的直线为渐近线.例2.1.观察数列}1{+n n 当∞→n 时的变化趋势.解:输入命令:n=1:100；xn=n./(n+1)得到该数列的前100项，从这前100项看出，随n 的增大，1+n n与1非常接近，画出nx 的图形.stem(n ，xn) 或for i=1:100;plot(n(i),xn(i),’r’) hold on end其中for … end 语句是循环语句,循环体内的语句被执行100次,n(i)表示n 的第i 个分量.由图可看出，随n 的增大，点列与直线1=y 无限接近，因此可得结论:1lim+∞→n nn =１.对函数的极限概念，也可用上述方法理解.例2.2.分析函数x x x f 1sin )(=，当0→x 时的变化趋势. 解:画出函数)(x f 在]1,1[-上的图形.x=-1:0.01:1；y=x.*sin(1./x)；plot(x ，y)从图上看，x x 1sin 随着x 的减小，振幅越来越小趋近于0，频率越来越高作无限次振荡.作出x y ±=的图象.hold on ；plot(x ，x ，x ，-x)例2.3.分析函数x x f 1sin )(=当0→x 时的变化趋势. 解:输入命令:x=-1:0.01:1；y=sin(1./x)；plot(x ，y)从图上看，当0→x 时，x 1sin 在-1和1之间无限次振荡，极限不存在.仔细观察该图象，发现图象的某些峰值不是１和-１，而我们知道正弦曲线的峰值是１和-１，这是由于自变量的数据点选取未必使x 1sin 取到１和-１的缘故，读者可试增加数据点，比较它们的结果． 例2.4.考察函数x xx f sin )(=当0→x 时的变化趋势. 解:输入命令:x=linspace(-2*pi ，2*pi ，100)；y=sin(x)./x ；plot(x ，y)从图上看，x xsin 在0=x 附近连续变化，其值与1无限接近，可见x xx sin lim 0→=1.例2.5.考察xx x f )11()(+=当∞→x 时的变化趋势. 解:输入命令:x=1:20:1000；y=(1+1./x).^x ；plot(x ，y)从图上看，当∞→x 时，函数值与某常数无限接近，我们知道，这个常数就是e . ３．求函数极限例2.6.求)1311(lim 31+-+-→x x x .解:输入命令:syms x ；f=1/(x+1)-3/(x^3+1)；limit(f ，x ，-1) 得结果ans=-1.画出函数图形.ezplot(f)；hold on ；plot(-1，-1，’r.’)例2.7.求30sin limx xtgx x -→.解:输入命令:limit((tan(x)-sin(x))/x^3) 得结果:ans=1/2例2.8.求xx x x )11(lim -+∞→.解:输入命令:limit(((x+1)/(x-1))^x ，inf) 得结果:ans=exp(2)例2.9.求xx x +→0lim . 解:输入命令:limit(x^x ，x ，0，’right’) 得结果:ans =1 例2.10.求x ctgx 10)(lim+→.解:输入命令:limit((cot(x))^(1/log(x))，x ，0，’right’) 得结果:ans=exp(-1)４．求方程的解. 例2.11.解方程02=++c bx ax . 解:输入命令:syms a b c x ；f=a*x^2+b*x+c ；solve(f) 得结果:ans=[ 1/2/a*(-b+(b^2-4*a*c)^(1/2))] [ 1/2/a*(-b-(b^2-4*a*c)^(1/2))]如果不指明自变量，系统默认为x ，也可指定自变量，比如指定b 为自变量. solve(f ，b)得结果:ans=-(a*x^2+c)/x例2.12.解方程0155=--x x . 解:输入命令:f=x^5-5*x-1；solve(f) 得结果:ans=[ -1.4405003973415600893186320629653] [ -.20006410262997539129073370075959] [.49456407933505360591791681140791e-1 -1.4994413672391491358223492788056*i] [.49456407933505360591791681140791e-1 +1.4994413672391491358223492788056*i] [ 1.5416516841045247594257824014433] 画出图象ezplot(f ，[-2，2])；hold on ；plot([-2，2]，[0，0])实验3 导数及偏导数计算实验目的１．进一步理解导数概念及其几何意义. ２．学习matlab 的求导命令与求导法.实验内容１．学习matlab 命令.建立符号变量命令sym 和syms 调用格式: x=sym(｀x ｀)， 建立符号变量x ；syms x y z ， 建立多个符号变量x ，y ，z ； matlab 求导命令diff 调用格式:diff (函数)(x f ) ， 求)(x f 的一阶导数)(x f '； diff (函数)(x f ， n ) ， 求)(x f 的n 阶导数)()(x fn (n 是具体整数)；diff (函数),(y x f ，变量名x )， 求),(y x f 对x 的偏导数x f∂∂；diff (函数),(y x f ， 变量名x ，n ) ，求),(y x f 对x 的n 阶偏导数n n x f∂∂；matlab 求雅可比矩阵命令jacobian ，调用格式:jacobian ([函数),,(z y x f ；函数),,(z y x g ； 函数),,(z y x h ]， [z y x ,,])给出矩阵:⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂z h yh xh z g y g x g z f y f x f２．导数概念.导数是函数的变化率，几何意义是曲线在一点处的切线斜率. （１）点导数是一个极限值.例3.1.设xe xf =)(，用定义计算)0(f '.解:)(x f 在某一点0x 的导数定义为极限:x x f x x f x ∆-∆+→∆)()(lim000我们记x h ∆=，输入命令:syms h ；limit((exp(0+h)-exp(0))/h ，h ，0) 得结果:ans=１．可知1)0(='f（２）导数的几何意义是曲线的切线斜率.例3.2.画出xe xf =)(在0=x 处()1,0(P )的切线及若干条割线，观察割线的变化趋势.解:在曲线x e y =上另取一点),(he h M ，则PM 的方程是:011--=--h e x y h .即 11+-=x h e y h取01.0,1.0,1,2,3=h ，分别作出几条割线.h=[3，2，1，0.1，0.01]；a=(exp(h)-1)./h ；x=-1:0.1:3； plot(x ，exp(x)，’r’)；hold on for i=1:5；plot(h(i)，exp(h(i))，’r.’) plot(x ，a(i)*x+1) endaxis square作出xe y =在0=x 处的切线1+=x yplot(x ，x+1，’r’)从图上看，随着M 与P 越来越接近，割线PM 越来越接近曲线的割线． ３．求一元函数的导数.（１）)(x f y =的一阶导数.例3.3.求x x y )sin(=的导数.解:打开matlab 指令窗，输入指令:dy_dx=diff(sin(x)/x). 得结果:dy_dx=cos(x)/x-sin(x)/x^2.matlab 的函数名允许使用字母、空格、下划线及数字，不允许使用其他字符，在这里我们用dy_dx 表示x y '. 例3.4.求)ln(sin x y =的导数.解: 输入命令:dy_dx=diff(log(sin(x))). 得结果:dy_dx=cos(x)/sin(x).在matlab 中，函数x ln 用log(x)表示，而log10(x)表示x lg . 例3.5.求202)2(x x y +=的导数.解: 输入命令:dy_dx=diff((x^2+2*x)^20). 得结果:dy_dx=20*(x^2+2*x)^19*(2*x+2).注意x 2输入时应为2*x.例3.6.求xx y =的导数. 解: 输入命令:dy_dx=diff(x^x). 得结果:dy_dx =x^x*(log(x)+1).利用matlab 命令diff 一次可以求出若干个函数的导数. 例3.7.求下列函数的导数:1.5221+-=x x y .2.x x y 2cos 2cos 22+=. 3.xy sin 34=.4.x y ln ln 4=. 解: 输入命令:a=diff([sqrt(x^2- 2*x+5)，cos(x^2)+2*cos(2*x)，4^(sin(x))， log(log(x))]). 得结果: a=[1/2/(x^2-2*x+5)^(1/2)*(2*x-2)，-2*sin(x^2)*x-4*sin(2*x)， 4^sin(x)*cos(x)*log(4)， 1/x/log(x)]. dy1_dx=a(1)↵.dy1_dx=1/2/(x^2-2*x+5)^(1/2)*(2*x-2). dy2_dx=a(2)↵.dy2_dx=-2*sin(x^2)*x-4*sin(2*x). dy3_dx=a(3)↵.dy3_dx=4^sin(x)*cos(x)*log(4).dy4_dx=a(4)↵.dy4_dx=1/x/log(x).由本例可以看出，matlab 函数是对矩阵或向量进行操作的，a(i)表示向量a 的第i 个分量.（２）参数方程所确定的函数的导数.设参数方程⎩⎨⎧==)()(t y y t x x 确定函数)(x f y =，则y 的导数)()(t x t y dx dy ''=. 例3.8.设⎩⎨⎧-=-=)cos 1()sin (t a y t t a x ，求dx dy .解: 输入命令:dx_dt=diff(a*(t-sin(t)))；dy_dt=diff(a*(1-cos(t)))； dy_dx=dy_dt/dx_dt. 得结果:dy_dx=sin(t)/(1-cos(t)). 其中分号的作用是不显示结果. ４．求多元函数的偏导数.例3.9.设 u=222z y x ++求 u 的一阶偏导数.解: 输入命令:diff((x^2+y^2+z^2)^(1/2)， x). 得结果:ans=1/(x^2+y^2+z^2)^(1/2)*x. 在命令中将末尾的x 换成y 将给出y 的偏导数:ans=1/(x^2+y^2+z^2)^(1/2)*y. 也可以输入命令:jacobian((x^2+y^2+z^2)^(1/2)，[x y]). 得结果:ans=[1/(x^2+y^2+z^2)^(1/2)*x ， 1/(x^2+y^2+z^2)^(1/2)*y]给出矩阵⎪⎭⎫ ⎝⎛∂∂∂∂y u x u ,.例3.10.求下列函数的偏导数:1.)(1x yarctg z =. 2.yx z =2.解: 输入命令:diff(atan(y/x). 得结果:ans=-y/x^2/(1+y^2/x^2). 输入命令:diff(atan(y/x)， y). 得结果:ans=1/x/(1+y^2/x^2). 输入命令:diff(x^y ， x). 得结果:ans=x^y*y/x . 输入命令:diff(x^y ， y). 得结果:ans=x^y*log(x).使用jacobian 命令求偏导数更为方便. 输入命令:jacobian([atan(y/x)，x^y]，[x ，y]). 得结果:ans=[ -y/x^2/(1+y^2/x^2)， 1/x/(1+y^2/x^2)] [ x^y*y/x ， x^y*log(x)]. ５．求高阶导数或高阶偏导数.例3.11.设xe x xf 22)(=，求)()20(x f.解:输入指令:diff(x^2*exp(2*x)，x ，20). 得结果: ans =99614720*exp(2*x)+20971520*x*exp(2*x)+1048576*x^2*exp(2*x)例3.12.设224623y x y x z +-=，求y x z y z x z ∂∂∂∂∂∂∂22222,,.解:输入命令:diff(x^6-3*y^4+2*x^2*y^2，x ，2)可得到22x z ∂∂:ans=30*x^4+4*y^2.将命令中最后一个x 换为y 得22y z ∂∂:ans=-36*y^2+4*x^2. 输入命令:diff(diff(x^6-3*y^4+2*x^2*y^2，x)，y)可得y x z ∂∂∂2:ans=8*x*y同学们可自己计算x y z∂∂∂2比较它们的结果.注意命令:diff(x^6-3*y^4+2*x^2*y^2，x ，y)，是对y 求偏导数，不是求y x z∂∂∂2.６．求隐函数所确定函数的导数或偏导数 例3.13.设e e x xy =+-ln ，求dx dy解:e ex y x F xy -+=-ln ),(，先求x F '，再求y F '. 输入命令:df_dx=diff(log(x)+exp(-y/x)-exp(1)，x)得到x F ':df_dx=1/x+y/x^2*exp(-y/x). 输入命令:df_dy=diff(log(x)+exp(-y/x)-exp(1)，y) 得到y F ':df_dy=-1/x*exp(-y/x) 输入命令:dy_dx=-df_dx/df_dy 可得所求结果:dy_dx=-(-1/x-y/x^2*exp(-y/x))*x/exp(-y/x).例3.14.设0)()cos()sin(=++xz tg yz xy ，求y zx z ∂∂∂∂,解:)()cos()sin()(xz tg yz xy x F ++=输入命令:a=jacobian(sin(x*y)+cos(y*z)+tan(z*x)，[x ，y ，z])可得矩阵()z y x F F F ''',,a=[cos(x*y)*y+(1+tan(z*x)^2)*z ，cos(x*y)*x-sin(y*z)*z ， -sin(y*z)*y+(1+tan(z*x)^2)*x]. 输入命令:dz_dx=-a(1)/a(3) 得:dz_dx=(-cos(x*y)*y-(1+tan(z*x)^2)*z)/(-sin(y*z)*y+(1+tan(z*x)^2)*x) 输入命令:dz_dy=-a(2)/a(3) 得:dz_dy=(-cos(x*y)*x+sin(y*z)*z)/(-sin(y*z)*y+(1+tan(z*x)^2)*x)实验4 积分计算实验目的１．通过本实验加深理解积分理论中分割、近似、求和、取极限的思想方法. ２．学习并掌握用matlab 求不定积分、定积分、二重积分、曲线积分的方法. ３．学习matlab 命令sum 、symsum 与int . 实验内容１．学习matlab 命令（１）求和命令sum 调用格式.sum(x)，给出向量x 的各个元素的累加和，如果x 是矩阵，则sum(x)是一个元素为x 的每列列和的行向量.例4.1.x=[1，2，3，4，5，6，7，8，9，10]；↵ sum(x)↵ ans=55例4.2.x=[1，2，3；4，5，6；7，8，9]↵ x=1 2 3 4 5 6 7 8 9 sum(x)↵ans=12 15 18（２）求和命令symsum 调用格式.symsum(s ，n)， 求∑nssymsum(s ，k ，m ，n)，求∑=nm k s当x 的元素很有规律，比如为表达式是)(k s 的数列时，可用symsum 求得x 的各项和，即 symsum ),1),((n k s =)()2()1(n s s s +++symsum )()1()(),,),((n s m s m s n m k k s ++++= 例4.3.syms k n ↵ symsum(k ，1，10)↵ ans=55symsum(k^2，k ，1，n)↵ans=1/3*(n+1)^3-1/2*(n+1)^2+1/6*n+1/6 （３）matlab 积分命令int 调用格式int (函数)(x f ) 计算不定积分⎰dx x f )(int (函数),(y x f ，变量名x ) 计算不定积分⎰dx y x f ),(int (函数b a x f ,),() 计算定积分⎰badxx f )(int (函数),,(y x f 变量名b a x ,,) 计算定积分⎰badxy x f ),(２．计算不定积分 例4.4.计算xdxx ln 2⎰解:输入命令:int(x^2*log(x)) 可得结果:ans=1/3*x^3*log(x)-1/9*x^3 注意设置符号变量.例4.5.计算下列不定积分: 1.dxx a ⎰-222.⎰++dx x x 31313.⎰xdxx arcsin 2解:首先建立函数向量. syms xsyms a realy=[sqrt(a^2-x^2)，(x-1)/(3*x-1)^(1/3)，x^2*asin(x)]； 然后对y 积分可得对y 的每个分量积分的结果.int(y ，x)↵ ans =[1/2*x*(a^2-x^2)^(1/2)+1/2*a^2*asin((1/a^2)^(1/2)*x)， -1/3*(3*x-1)^(2/3)+1/15*(3*x-1)^(5/3)，1/3*x^3*asin(x)+1/9*x^2*(1-x^2)^(1/2)+2/9*(1-x^2)^(1/2)] ３．定积分的概念.定积分是一个和的极限.取xe xf =)(，积分区间为]1,0[，等距划分为20个子区间. x=linspace(0，1，21)；选取每个子区间的端点，并计算端点处的函数值. y=exp(x)；取区间的左端点乘以区间长度全部加起来. y1=y(1:20)；s1=sum(y1)/20 s1=1.6757s1可作为⎰10dx e x 的近似值.若选取右端点:y2=y(2:21)；s2=sum(y2)/20 s2=1.7616s2也可以作为⎰10dx e x 的近似值.下面我们画出图象.plot(x ，y)；hold onfor i=1:20 fill([x(i)，x(i+1)，x(i+1)，x(i)，x(i)]，[0，0，y(i)，y(i)，0]，’b’)end如果选取右端点，则可画出图象. for i=1:20； fill([x(i)，x(i+1)，x(i+1)，x(i)，x(i)]，[0，0，y(i+1)，y(i+1)，0]，’b’) hold on endplot(x ，y ，’r’)在上边的语句中，for … end 是循环语句，执行语句体内的命令20次，fill 命令可以填充多边形，在本例中，用的是兰色(blue)填充.从图上看211s dx e s x <<⎰，当分点逐渐增多时，12s s -的值越来越小，读者可试取50个子区间看一看结果怎样.下面按等分区间计算∑∑=∞→=∞→=∆ξn i n i n ni i n n ex f 111lim)(limsyms k ns=symsum(exp(k/n)/n ，k ，1，n)； limit(s ，n ，inf) 得结果ans=exp(1)-1４．计算定积分和广义积分.例4.6.计算⎰10dx e x .解:输入命令:int(exp(x)，0，1)得结果ans=exp(1)-1.这与我们上面的运算结果是一致的. 例4.7.计算⎰-201dxx解:输入命令:int(abs(x-1)，0，2)得结果ans =1.本例用mathematica 软件不能直接求解. 例4.8.判别广义积分⎰+∞11dx x p 、⎰+∞∞--πdx e x 2221与⎰-202)1(1dx x 的敛散性，收敛时计算积分值.解:对第一个积分输入命令:syms p real ；int(1/x^p ，x ，1，inf)得结果ans =limit(-1/(p-1)*x^(-p+1)+1/(p-1)，x = inf).由结果看出当1<p 时，x^(-p+1)为无穷，当1>p 时，ans=1/(p-1)，这与课本例题是一致的.对第二个积分输入命令:int(1/(2*pi)^(1/2)*exp(-x^2/2)，-inf ，inf) 得结果:ans=7186705221432913/18014398509481984*2^(1/2)*pi^(1/2) 由输出结果看出这两个积分收敛.对后一个积分输入命令:int(1/(1-x)^2，0，2)结果得ans=inf .说明这个积分是无穷大不收敛. 例4.9.求积分⎰t dxxx 0sin解:输入命令:int(sin(x)/x ，0，t)，可得结果sinint(t)，通过查帮助(help sinint )可知sinint(t)=⎰t dxxx 0sin，结果相当于没求!实际上matlab 求出的只是形式上的结果，因为这类积分无法用初等函数或其值来表示.尽管如此，我们可以得到该函数的函数值.输入vpa(sinint(0.5))可得sinint(0.5)的值.５．二重积分计算 例4.10.求二次积分⎰⎰+12102x x xydydx解:输入命令:int(int(x*y ，y ，2*x ，x^2+1)，x ，0，1) 得结果ans=1/12. 例4.11.求⎰⎰≤++π12222))(sin(y x dxdyy x解:积分区域用不等式可以表示成2211,11x y x x -≤≤--≤≤-，二重积分可化为二次积分⎰⎰----+π22112211)(sin(x x dyy x dx，输入命令:int(int(sin(pi*(x^2+y^2))，y ，-sqrt(1-x^2)，sqrt(1-x^2))，x ，-1，1) 由输出结果可以看出，结果中仍带有int ，表明matlab 求不出这一积分的值.采用极坐标可化为二次积分⎰⎰ππ20102)sin(drr r da ，输入命令:int(int(r*sin(pi*r^2)，r ，0，1)，a ，0，2*pi) 可得结果为ans=2.６．曲线积分例4.12.求曲线积分⎰L xyds ，其中L 为曲线122=+y x 在第一象限内的一段.解:曲线的参数方程是)20(,sin ,cos π≤≤==t t y t x 曲线积分可以化为⎰π⋅0s i n c o s t d tt .输入命令:int(cos(t)*sin(t)，0，pi/2) 执行后即可求出曲线积分结果1/2.实验5 matlab 自定义函数与导数应用实验目的１．学习matlab 自定义函数.２．加深理解罗必塔法则、极值、最值、单调性. 实验内容１．学习matlab 自定义函数及求函数最小值命令.函数关系是指变量之间的对应法则，这种对应法则需要我们告诉计算机，这样，当我们输入自变量时，计算机才会给出函数值，matlab 软件包含了大量的函数，比如常用的正弦、余弦函数等.matlab 允许用户自定义函数，即允许用户将自己的新函数加到已存在的matlab 函数库中，显然这为matlab 提供了扩展的功能，无庸置疑，这也正是matlab 的精髓所在.因为matlab 的强大功能就源于这种为解决用户特殊问题的需要而创建新函数的能力.matlab 自定义函数是一个指令集合，第一行必须以单词function 作为引导词，存为具有扩展名“.m ”的文件，故称之为函数M －文件.函数M －文件的定义格式为：function 输出参数＝函数名（输入参数） 函数体 …… 函数体一旦函数被定义，就必须将其存为M －文件，以便今后可随时调用.比如我们希望建立函数12)(2++=x x x f ，在matlab 工作区中输入命令：syms x ；y=x^2+2*x+1；不能建立函数关系，只建立了一个变量名为y 的符号表达式，当我们调用y 时，将返回这一表达式.y ↵y=x^2+2*x+1当给出x 的值时，matlab 不能给出相应的函数值来. x=3；y ↵y=x^2+2*x+1 如果我们先给x 赋值.x=3；y=x^2+2*x+1 得结果：y=16若希望得出2|=x y 的值，输入： x=2；y ↵得结果：y=16，不是2=x 时的值.读者从这里已经领悟到在matlab 工作区中输入命令：y=x^2+2*x+1不能建立函数关系，如何建立函数关系呢？我们可以点选菜单Fill\New\M-fill 打开matlab 文本编辑器，输入： function y=f1(x) y=x^2+2*x+1；存为f1.m .调用该函数时，输入： syms x ；y=f1(x)↵ 得结果：y= x^2+2*x+1.输入: y1=f1(3)↵ 得结果：y1=16matlab 求最小值命令fmin 调用格式：fmin(‘fun’，a ，b) 给出)(x f 在),(b a 上的最小值点. ２．自定义函数例5.1.建立正态分布的密度函数22)(21),.,(μ--σπ=μσx e x f解：打开文本编辑器，输入：function y=zhengtai(x ，a ，b)y=1/sqrt(2*pi)/a*exp(-(x-b).^2/2/a^2)； 存为zhengtai.m .调用时可输入命令： y=zhengtai(1，1，0)得结果：y=0.2420.此即)0,1,1(f 的值.如果想画出标准正态分布的密度函数的图象，输入： ezplot(zhengtai(x ，1，0))例5.2.解一元二次方程02=++c bx ax .解：我们希望当输入c b a ,,的值时，计算机能给出方程的两个根.在文本编辑器中建立名为rootquad.m 的文件.function [x1，x2]=rootquad(a ，b ，c) d=b*b-4*a*c ；x1=(-b+sqrt(d))/(2*a) x2=(-b-sqrt(d))/(2*a) 比如求方程07322=-+x x 的根，可用语句： [r1，r2]=rootquad(2，3，-7)得结果：r1=1.2656 r2=-2.7656２．验证罗必塔法则.罗必塔法则是指在求00及∞∞的极限时，可用导数之比的极限来计算（如果导数之比的极限存在或∞）例5.3.以x ba xx x -→0lim 为例验证罗必塔法则.解：这是00型极限f=a^x-b^x ；g=x ；L=limit(f/g ，x ，0)得结果：L=log(a)-log(b)df=diff(f ，x)；dg=diff(g ，x)；L1=limit(df/dg ，x ，0) 得结果：L1=log(a)-log(b) 从结果看出：L=L1，即x b a x x x -→0lim=x b a x x x '-→'0)(lim４．函数的单调性与极值.例5.4.求函数396)(23++-=x x x x f 的单调区间与极值.解：求可导函数的单调区间与极值，就是求导函数的正负区间与正负区间的分界点，利用matlab 解决该问题，我们可以先求出导函数的零点，再画出函数图象，根据图象可以直观看出函数的单调区间与极值．输入命令：f=x^3-6*x^2+9*x+3；df=diff(f ，x)；s=solve(df) 得结果：ans=[1，3]，画出函数图象． ezplot(f ，[0，4])从图上看，)(x f 的单调增区间为)1,(-∞、),1(+∞，单调减区间是)3,1(，极大值7)1(=f ，极小值3)3(=f .我们可以建立一个名为dandiao.m 的M —文件，用来求求函数的单调区间. disp(‘输入函数（自变量为x ）’) syms xf=input('函数f(x)=') df=diff(f)； s=solve(df) a=[]；for i=1:size(s)； a(i)=s(i)； endezplot(f ，[min(a)-1，max(a)+1])要求函数)1ln(x x y +-=的单调区间与极值，可调用dandiao.m .输入： dandiao ↵在matlab 工作区出现以下提示： 输入函数（自变量为x ） 函数f(x)=在光标处输入：x-log(1+x)，可得结果s=0.从图上看，)(x f 的单调增区间为),0(+∞，单调减区间是)0,(-∞，极小值0)0(=f . ５．函数的最值调用求函数最小值命令fmin 时，可得出函数的最小值点，为求最小值，必须建立函数M —文件.例5.5.求函数1)3()(2--=x x f 在区间)5,0(上的最小值. 解：我们可以建立一个名为f.m 的函数M －文件. function y=f(x) y=(x-3).^2-1； 并且调用fminx=fmin((‘f’，0，5)得：x=3，)(x f 在最小值点处的值（函数最小值）是1)3(-=f .求最大值时可用x=fmin(‘-f(x)’，a ，b)实验6 matab 矩阵运算与数组运算实验目的：１．理解矩阵及数组概念．２．掌握matlab 对矩阵及数组的操作命令．实验内容：１．矩阵与数组的输入． 对于较小较简单的矩阵，从键盘上直接输入矩阵是最常用的数值矩阵创建方法．用这种方法输入矩阵时注意以下三点：（1）整个输入矩阵以方括号“[ ]”为其首尾； （2）矩阵的元素必须以逗号“，”或空格分隔； （3）矩阵的行与行之间必须用分号“；”或回车键隔离． 例１：下面的指令可以建立一个3行4列的矩阵 a ． a=[1 2 3 4;5 6 7 8;9 10 11 12] (下面是屏幕的显示结果) a =1 2 3 4 5 6 7 8 9 10 11 12 分号“；”有三个作用：（1）在“[ ]”方括号内时它是矩阵行间的分隔符．例子如上． （2）它可用作指令与指令间的分隔符．（3）当它存在于赋值指令后时，该指令执行后的结果将不显示在屏幕上． 例如，输入指令：b=[1 2 0 0;0 1 0 0;1 1 1 1];矩阵b 将不被显示，但b 已存放在matlab 的工作内存中，可随时被以后的指令所调用或显示．例如，输入指令： b得结果： b =1 2 0 0 0 1 0 0 1 1 1 1数值矩阵的创建还可由其他方法实现．如：利用matlab 函数和语句创建数值矩阵；利用m 文件创建数值矩阵；从其他文件获取数值矩阵．有兴趣的读者可参阅其他参考书．数组可以看成特殊的矩阵，即1行n列的矩阵，数组的输入可以采用上面矩阵的输入方法．例２：输入以下指令以建立数组c．c=[1 2 3 4 5 6 7 8]c =1 2 3 4 5 6 7 8另外还有两种方法输入数组．请看下面两个例子．例３：在0和2中间每隔0.1一个数据建立数组d．解：输入指令：d=0:0.1:2d =Columns 1 through 70 0.1000 0.2000 0.3000 0.4000 0.50000.6000Columns 8 through 140.7000 0.8000 0.9000 1.0000 1.1000 1.20001.3000Columns 15 through 211.4000 1.5000 1.6000 1.7000 1.8000 1.90002.0000注意“：”的使用方法.例４：在0和2之间等分地插入一些分点，建立具有10个数据点的数组e．解：输入指令：e=linspace(0,2,10)e =Columns 1 through 70 0.2222 0.4444 0.6667 0.8889 1.11111.3333Columns 8 through 101.5556 1.77782.0000linspace(a,b,n)将建立从a到b有n个数据点的数组.２．常用矩阵的生成．matlab为方便编程和运算，提供了一些常用矩阵的生成指令：n⨯单位矩阵eye(n) nn⨯全1矩阵ones(n) nn⨯零矩阵zeros(n) nm⨯标准型矩阵eye(m,n) nm⨯全1矩阵ones(m,n) nm⨯零矩阵zeros(m,n) neye(size(A)) 与A同型的标准型矩阵ones(size(A)) 与A同型的全1矩阵zeros(size(A))与A同型的零矩阵其中指令size(A)给出矩阵A的行数和列数．例５：生成以下矩阵．3⨯零矩阵．(1)33 全1矩阵．(2)6(3)与例1中矩阵a同型的标准型矩阵．解：输入下面指令：d=zeros(3)d =0 0 00 0 00 0 0e=ones(3,6)e =1 1 1 1 1 11 1 1 1 1 11 1 1 1 1 1f=eye(size(a))f =1 0 0 00 1 0 00 0 1 0３．矩阵元素的标识矩阵的元素、子矩阵可以通过标量、向量、冒号的标识来援引和赋值．(1)矩阵元素的标识方式A(ni,nj)．ni,nj都是标量．若它们不是整数，则在调用格式中会自动圆整到最临近整数．ni指定元素的行位置，nj指定元素的列位置．(2)子矩阵的序号向量标识方式A(v,w)．v,w是向量，v,w中的任意一个可以是冒号“：”，他表示取全部行(在v位置)或全部列(在w位置)．v,w中所用序号必须大于等于1且小于等于矩阵的行列数．例６：元素和矩阵的标识a=[1 2 3 4;5 6 7 8;9 10 11 12]a =1 2 3 45 6 7 89 10 11 12a24=a(2,4)a24 =8a1=a([1,2],[2,3,4])a1 =2 3 46 7 8a2=a([1,2],[2,3,1])a2 =2 3 16 7 5a3=a([3,1],:)a3 =9 10 11 121 2 3 4a([1,3],[2,4])=zeros(2)a =1 0 3 05 6 7 89 0 11 0４．矩阵运算和数组运算．矩阵运算的指令和意义如下：A' 矩阵A的共轭转置矩阵，当A是实矩阵时，A'是A的转置矩阵．A+B 两个同型矩阵A与B相加．A-B 两个同型矩阵A与B相减．A*B 矩阵A与矩阵B相乘，要求A的列数等于B的行数．s+B 标量和矩阵相加（matlab约定的特殊运算，等于s加B的每一个分量）．s-B B-s 标量和矩阵相减（matlab约定的特殊运算，含意同上）．s*A 数与矩阵A相乘．例７：a=[1 2 3;4 5 6]a =1 2 34 5 6b=[-1 0 1;3 1 2]b =-1 0 13 1 2a'ans =1 42 53 6a+bans =0 2 47 6 8a-bans =2 2 21 4 41+aans =2 3 45 6 7a-1ans =0 1 23 4 52*bans =-2 0 26 2 4c=[2 4;1 3;0 1]c =2 41 30 1a*cans =4 1313 37数组可以看成特殊矩阵即一行n列的矩阵，矩阵运算的指令和含意同样适用于数组运算．如果在运算符前加“．”，含意将有所不同．A.*B 同维数组或同型矩阵对应元素相乘．A./B A的元素被B的元素对应除．A.^n A的每个元素n次方．p.^A 以p为底，分别以A的元素为指数求幂．例８：a=[1 2 3;4 5 6]a =1 2 34 5 6b=[-1 0 1;3 1 2]b =-1 0 13 1 2a.*bans =-1 0 312 5 12a./bWarning: Divide by zero.ans =-1.0000 Inf 3.00001.3333 5.0000 3.0000a.^2ans =1 4 916 25 362.^aans =2 4 816 32 64实验7矩阵与线性方程组实验目的：１．掌握matlab求矩阵的秩命令．２．掌握matlab求方阵的行列式命令．３．理解逆矩阵概念，掌握matlab求逆矩阵命令．４．会用matlab求解线性方程组．实验内容：１．矩阵的秩.指令rank(A)将给出矩阵A的秩．例１：a=[3 2 -1 -3 -2;2 -1 3 1 -3;7 0 5 -1 -8]a =3 2 -1 -3 -22 -13 1 -37 0 5 -1 -8rank(a)ans =2２．方阵的行列式.指令det(A)给出方阵A的行列式.例２：b=[1 2 3 4;2 3 4 1;3 4 1 2;4 1 2 3];det(b)ans =160det(b')ans =160c=b;c(:,1)=2*b(:,1);det(c)ans =320det(b(:,[3 2 1 4]))ans =-160d=b;d(2,:);det(d)ans =160你能解释上例中的运算结果吗?在这里我们实际上验证了行列式的性质．３．逆矩阵指令inv(A)给出方阵A的逆矩阵，如果A不可逆，则inv(A)给出的矩阵的元素都是Inf．例３：设⎪⎪⎪⎭⎫⎝⎛=343122321A，求A的逆矩阵．解：输入指令：A=[1 2 3;2 2 1;3 4 3];B=inv(A)B =1.0000 3.0000 -2.0000-1.5000 -3.0000 2.50001.0000 1.0000 -1.0000还可以用伴随矩阵求逆矩阵，打开m 文件编辑器，建立一个名为companm 的M-文件文件内容为：function y=companm(x)[n,m]=size(x);y=[];for j=1:n;a=[];for i=1:n;x1=det(x([1:i-1,i+1:n],[1:j-1,j+1:n]))*(-1)^(i+j);a=[a,x1];endy=[y;a];end利用该函数可以求出一个矩阵的伴随矩阵．输入命令：C=1/det(A)*companm(A)C =1.0000 3.0000 -2.0000-1.5000 -3.0000 2.50001.0000 1.0000 -1.0000利用初等变换也可以求出逆矩阵，构造n 行2n 列的矩阵(A E)，并进行行初等变换，当把A 变为单位矩阵时，E 就变成了A 的逆矩阵．利用matlab 命令rref 可以求出矩阵的行简化阶梯形．输入命令：D=[A,eye(3)]D =1 2 3 1 0 02 2 1 0 1 03 4 3 0 0 1rref(D)ans =1.0000 0 0 1.0000 3.0000 -2.0000 0 1.0000 0 -1.5000 -3.0000 2.5000 0 0 1.0000 1.0000 1.0000 -1.0000n m ⨯线性方程组B AX =的求解是用矩阵除来完成的,B A X \=，当n m =且A 可逆时，给出唯一解．这时矩阵除B A \相当于B A inv *)(；当m n >时，矩阵除给出方程的最小二乘解；当m n <时，矩阵除给出方程的最小范数解．例４：解方程组：⎪⎪⎩⎪⎪⎨⎧=-+=++=+-+=++-12121243132143214321x x x x x x x x x x x x x x 解：输入命令：a=[1 -1 1 2;1 1 -2 1;1 1 1 0;1 0 1 -1];b=[1;1;2;1];x=a\bx =0.83330.75000.41670.2500输入命令：z=inv(a)*bz =0.83330.75000.41670.2500例５：解方程组：⎪⎩⎪⎨⎧=-++-=-++-=--++8343242222543215432154321x x x x x x x x x x x x x x x解：方程的个数和未知数不相等，用消去法，将增广矩阵化为行简化阶梯形，如果系数矩阵的秩不等于增广矩阵的秩，则方程组无解；如果系数矩阵的秩等于增广矩阵的秩，则方程组有解，方程组的解就是行简化阶梯形所对应的方程组的解．输入命令：a=[2 1 1 -1 -2 2;1 -1 2 1 -1 4;2 -3 4 3 -1 8];rref(a)ans =1 0 0 0 0 00 1 0 -1 -1 00 0 1 0 -1 2由结果看出，4x ,5x 为自由未知量，方程组的解为：01=x542x x x +=532x x +=例６：解方程组：⎪⎪⎩⎪⎪⎨⎧=+--=--=-+-=+--0320030432142143214321x x x x x x x x x x x x x x x解：输入命令：a=[1 -1 -1 1;1 -1 1 -3;1 -1 0 -1;1 -1 -2 3];rref(a)ans =1 -1 0 -10 0 1 -20 0 0 00 0 0 0由结果看出，2x ,4x 为自由未知量，方程组的解为：421x x x +=432x x =实验8 常微分方程与级数实验目的：１．学习用matlab求解微分方程命令dsolve.２．学习matlab泰勒级数展开命令.３．巩固幂级数的收敛半径、和等概念.实验内容：１．学习matlab命令.matlab求解微分方程命令dsolve，调用格式为：dsolve(‘微分方程’)给出微分方程的解析解，表示为t的函数.dsolve(‘微分方程’，‘初始条件’)给出微分方程初值问题的解，表示为t的函数.dsolve(‘微分方程’，‘变量x’)给出微分方程的解析解，表示为x的函数.dsolve(‘微分方程’，‘初始条件’，‘变量x’)给出微分方程初值问题的解，表示为x的函数.求已知函数的taylor展开式taylor命令，调用格式为：taylor(函数f(x)) f(x)的5次taylor多项式.taylor(函数f(x)，n) f(x)的n-1次taylor多项式.taylor(函数f(x)，a) f(x)在a点的taylor多项式.求级数的和命令symsum调用格式为：symsum(S,n)，求∑n Ssymsum(S,k,m,n),求∑=nmkSmatlab求极限命令limit调用格式为：limit(函数f(x)，变量x，自变量的趋向值)２．求解一阶微分方程.微分方程在输入时，y'应输入Dy，y''应输入D2y等，D应大写．例１：求微分方程22xxexydxdy-=+的通解．解：输入命令：dsolve('Dy+2*x*y=x*exp(-x^2)')结果为ans =1/2*(1+2*exp(-2*x*t)*C1*exp(x^2))/exp(x^2)系统默认的自变量是t，显然系统把x当作常数，把y当作t的函数求解．输入命令：dsolve('Dy+2*x*y=x*exp(-x^2)','x')得正确结果：ans =1/2*(x^2+2*C1)/exp(x^2)例２：求微分方程0=-+'x e y y x 在初始条件e y x 21==下的特解．解：输入命令：dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')得结果为：ans =1/x*(exp(x)+exp(1))例３：求微分方程0cos 2)1(2=-+-x xy dx dy x 在初始条件10==x y 下的特解．解：输入命令：dsolve('(x^2-1)*Dy+2*x*y-cos(x)=0','y(0)=1','x')得结果为ans =1/(x^2-1)*(sin(x)-1)３．求解二阶微分方程．例４：求03=+'+''x e y y 的通解．解：输入命令：dsolve('D2y+3*Dy+exp(x)=0','x')得结果：ans =-1/4*exp(x)+C1+C2*exp(-3*x)例５：求解微分方程．02='-''y e y y 解：输入命令：dsolve('D2y-exp(2*y)*Dy=0','x')得结果：ans =1/2*log(-2*C1/(-1+exp(2*x*C1+2*C2*C1)))+x*C1+C2*C1４．taylor 展开式．例６：求函数y=cosx 在x=0点处的5阶taylor 展开式及在3π=x 处的6阶taylor 展开式．解：输入命令：syms x;taylor(cos(x))得结果：ans =1-1/2*x^2+1/24*x^4输入命令：taylor(cos(x),pi/3,7)得结果：ans =1/2-1/2*3^(1/2)*(x-1/3*pi)-1/4*(x-1/3*pi)^2+1/12*3^(1/2)*(x- 1/3*pi)^3+1/48*(x-1/3*pi)^4-1/240*3^(1/2)*(x-1/3*pi)^5-1/1440*(x-1/3*pi)^6 ５．级数求和．例７：求∑∞=121n n ．解：输入命令：syms n;symsum(1/2^n,1,inf) 得结果：ans =1例８：求幂级数∑∞=⨯12nnnnx的和函数．解：输入命令：symsum(x^n/(n*2^n),n,1,inf) 得结果ans =-log(1-1/2*x)例９：求幂级数∑∞=1nnnx的和函数．解：输入命令：symsum(n*x^n,n,1,inf) 得结果：ans =x/(x-1)^2６．判别级数敛散性．例１０：判断数项级数∑∞=+1)1(1nnn的收敛性．解：输入求和命令：symsum(1/(n*(n+1)),n,1,inf) 得结果：ans =1求和得是1，说明该级数收敛．例１１：判别级数∑∞=+1)1(sinnnnπ的敛散性．解：输入命令：symsum(sin(pi/(n*(n+1))),1,inf)得结果：ans =sum(sin(pi/n/(n+1)),pi = 1 .. inf)由执行结果看出仍含有sum，说明用matlab不能求出其和，可采用比较判别法，取比较级数为P级数∑∞=121nn,取二者通项比值的极限．输入命令：limit(sin(pi/(n*(n+1)))/(1/n^2),n,inf)得结果：ans =pi得值为pi,由所取P级数收敛，得知所要判别的级数也收敛．例１２：判别级数∑∞=⎪⎭⎫⎝⎛143nnn的敛散性．解：用比值判别法，输入命令：limit((n+1)*(3/4)^(n+1)/(n*(3/4)^n),n,inf)得结果：ans =3/4极限值小于1，由比值判别法知级数收敛．实际上输入求和命令：symsum(n*(3/4)^n,n,1,inf)得结果：ans =12。
数学实验报告的总结
数学实验报告的总结引言本次数学实验旨在验证金典的两数相加的方法是否真的可以得到正确的结果。
通过实际操作，我们将对该方法进行验证，并分析其准确性和可靠性。
实验步骤1. 构造实验样本：我们选择了一组随机数作为实验样本，以确保实验结果的客观性。
2. 执行两数相加操作：我们按照金典的两数相加方法进行运算。
3. 记录实验结果：将实验数据记录在表格中，以便后续对比和分析。
实验结果经过实验操作和数据统计，我们得到了以下结果：数字1 数字2 预测结果实际结果- -1 2 3 34 7 11 119 5 14 143 8 11 116 2 8 8结果分析通过对比实验结果，我们可以发现，预测结果和实际结果完全一致，证明金典的两数相加方法确实是可行的。
在本次实验中，我们使用的样本数较小，但结果的一致性和重复性使我们对该方法产生了足够的信心。
金典的两数相加方法是一种简洁而有效的方法，它可以在不需要大量计算和复杂操作的情况下，快速得到准确的结果。
这使得它非常适用于计算机程序设计中，特别是在需要频繁进行相加操作的场景下。
结论通过本次实验，我们验证了金典的两数相加方法的准确性和可靠性。
该方法在实际操作中得到了正确的结果，证明了它的可靠性。
我们可以将该方法应用于实际计算中，以便提高计算效率和准确性。
然而，需要注意的是，该方法只适用于较小的数字相加，当数字较大时可能会产生溢出的问题。
因此，在实际应用中，我们需要根据具体场景选择合适的相加方法，以确保计算的准确性和可靠性。
参考文献1. 金典. (2010). 数学计算方法. 清华大学出版社.2. 张三, 李四. (2015). 实验报告撰写指南. 中国教育出版社.致谢在此，我要对实验中帮助过我的同学和老师表示衷心的感谢。
没有你们的支持和帮助，本次实验的顺利进行和成功完成是不可能的。
谢谢大家的辛勤付出和支持！。
数值分析实验报告_清华大学_非线性方程的解法
非线性方程的解法实验1.算法设计与比较问题提出：非线性方程组的求解方法很多，基本的思想是线性化。
不同的方法效果如何，要靠计算的实践来分析、比较。
实验内容：考虑算法(1)牛顿法(2)拟牛顿法分别编写它们的matlab程序。
实验要求：(1)用上述方法，分别计算两个例子。
在达到精度相同的前提下，比较迭代次数、浮点运算次数和CPU时间等。
1.1程序清单为使用flops统计浮点运算次数，使用MATLAB5.3版本%f1.m原函数f1function y=f(x)y(1)=12*x(1)-x(2)^2-4*x(3)-7;y(2)=x(1)^2+10*x(2)-x(3)-8;y(3)=x(2)^3+10*x(3)-8;end%ff1.m原函数f1的雅克比矩阵function y=ff(x)y(1,:)=[12,-2*x(2),-4];y(2,:)=[2*x(1),10,-1];y(3,:)=[0,3*x(2)^2,10];end%f1.m原函数f2function y=f2(x)y(1)=3*x(1)-cos(x(2)*x(3)) -1/2;y(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;y(3)=exp(-x(1)*x(2))+20*x(3)+1/3*(10*pi-3);end%ff2.m原函数f2的雅克比矩阵function y=ff2(x)y(1,:)=[3,x(3)*sin(x(2)*x(3)),x(2)*sin(x(2)*x(3))];y(2,:)=[2*x(1),-2*81*(x(2)+0.1),cos(x(3))];y(3,:)=[-x(2)*exp(-x(1)*x(2)),-x(1)*exp(-x(1)*x(2)),20]; end%牛顿法(以第一个方程组为例)clear;x0=[0,0,0]';n=10;tol=1e-6;x(:,1)=x0;i=1;u=[1,1,1]';tic;while (norm(u)>tol*norm(x(:,i))&(i<n))A=ff1(x(:,i));b=f1(x(:,i))';u=-A\b;x(:,i+1)=x(:,i)+u;i=i+1;end;x(:,i)iter=i-1t=toc%拟牛顿法(以第一个方程组为例)clear;x0=[0,0,0]';n=10;tol=1e-6;x(:,1)=x0;i=1;p=[1,1,1]';A=ff1(x(:,1));tic;while (norm(p)>tol*norm(x(:,i))&(i<n))x(:,i+1)=x(:,i)-A\f1(x(:,i))';p=x(:,i+1)-x(:,i);q=f1(x(:,i+1))'-f1(x(:,i))';A=A+(q-A*p)*p'/norm(p,2)^2;i=i+1;end;iter=i-1t=tocx(:,i)1.2运行结果1.2.1第一个方程组精确解为*T =(0.886020214719037, 0.796444775323146, 0.749479574122230)x 取最大迭代次数n=5000，相对误差限Tol=1e-6 (1)取()(0)1,1,1x T=牛顿迭代法迭代3次收敛，浮点运算次数为440，每次迭代平均浮点运算次数为147，CPU 耗时t =0(s)拟牛顿法迭代4次收敛，浮点运算次数为1048，每次迭代平均浮点运算次数为262，CPU 耗时t =0(s)(2)取()(0)000x T =，， 牛顿迭代法迭代4次收敛，浮点运算次数为510，每次迭代平均浮点运算次数为128，CPU 耗时t =1.600e-002(s)拟牛顿法迭代6次收敛，浮点运算次数为1493，每次迭代平均浮点运算次数为248，CPU 耗时t =1.50e-002(s)(3)取()(0)50,5050x T=，牛顿迭代法迭代15次收敛，浮点运算次数为2118，每次迭代平均浮点运算次数为141，CPU 耗时t =1.600e-002(s)拟牛顿法迭代338次收敛，浮点运算次数为88454，每次迭代平均浮点运算次数为262，CPU 耗时t =3.100e-002(s)1.2.2第二个方程组精确解为*T =(0.886020214719037, 0.796444775323146, 0.749479574122230)x 取最大迭代次数n=5000，相对误差限Tol=1e-6(1)取()(0)000x T=，， 牛顿迭代法迭代5次收敛，浮点运算次数为776，每次迭代平均浮点运算次数为155.2，CPU 耗时t =0(s)拟牛顿法迭代6次收敛，浮点运算次数为1635，每次迭代平均浮点运算次数为273，CPU 耗时t =0(s)(2)取()(0)888x T=，， 牛顿迭代法迭代9次收敛，浮点运算次数为1519，每次迭代平均浮点运算次数为169，CPU 耗时t =0(s)拟牛顿法迭代21次收敛，浮点运算次数为5924，每次迭代平均浮点运算次数为282，CPU 耗时t =1.600e-002(s)(3)对于离精确解更远的初值（如()(0)101010x T=，，），在计算中会出现奇异或接近奇异的矩阵，计算结果误差很大或计算根本无法进行下去。
清华杨顶辉数学实验六作业
Y= 1.0e-006 *
-0.0001
0.0133
0.0780
-0.1523
>> XT0=[0,0.9,0,70]; [XT,Y]=fsolve(@azeofun,XT0,[],n,P,a,b,c,Q) XT = 0.0000 0.7803 0.0000 76.9613
Y= 1.0e-008 * -0.0001 0.0446 结论：列表如下： XT0 [1/3,1/3,1/3,50] [0,0.9,0,70] x1 0 0 -0.1800 x2 0.5858 0.7803 -0.0019 x3 0.4142 0 x4 0 0.2197 T 71.9657 76.9613
6-8，假设商品在 t 时刻的市场价格为 p（t） ，需求函数为 D(p(t))=c-dp(t)（c，d≥0） 。而生产 方的期望价格为 q（t） ，供应函数为 S（q（t） ） 。当供销平衡时 S(q(t))=D（q(t)） 。若期望价 格与市场价格不符，商品市场不均衡，生产方 t+1 时期的期望价格将会调整，方式为 q(t+1)-q(t)=r[p(t)-q(t)](0<r<1)，以 p(t)=[c-D(p(t))]/d=[c-S(q(t))]/d 代入，得到关于 q(t)的递推方 程。设 S(x)=arctan(μx)，μ=4.8，d=0.25，r=0.3，以 c 为可变参数，讨论期望价格 q(t)的变 化规律，是否有混沌现象产生？并找出几个分叉点，观察分叉点极限趋势是否符合 Feigenbaum 常数揭示的规律。 解：模型建立：由题目可知： q(t+1)=rc/d-(arctan μq(t))r/d+(1-r)q(t)，代入数据得： q(t+1)=1.2c-1.2arctan4.8q(t)+0.7q(t) function y=iter01(x,c) y=1.2*c-1.2*arctan(4.8*x)+0.7*x; function chaos(iter_fun,x0,r,n) kr=0; for rr=r(1):r(3):r(2) kr=kr+1; y(kr,1)=feval(iter_fun,x0,rr); for i=2:n(2) y(kr,i)=feval(iter_fun,y(kr,i-1),rr); end end plot([r(1):r(3):r(2)],y(:,n(1)+1:n(2)),'k.');
清华大学分光计实验报告
'叉丝中心线平分，平行光管
○4 三棱镜两光学面法线垂直于分光计主轴 先粗调小平台使之平行于度盘平面 再以三棱镜一条边垂直两个螺钉连线的方位放置三棱镜于小平台，通过
望远镜对准其中一个光学面观察 2-8 的像，调节小平台螺钉使之回到 2-8 位 置，在旋转小平台，对准另一个光学面重复上述操作，如此重复，直到两个 光学面都使像到达 2-8 位置。
2
3
游标Ⅰ 游标Ⅱ 游标Ⅰ 游标Ⅱ
519’ 18518’ 1941’ 19940’
24524’ 6522’ 25945’ 7943’
i T1 T2 11955’ 11955’ 11955’ 11956’ 11956’ 11957’
(I II ) / 2
11955’
11956’
11957’
计算 A
主要配件：三棱镜，平面镜
分光计外形如图 1 （2）分光计的调节
先进行粗调，通过自己的感觉，调整使得望远镜，小平台，平行光管与度盘 平面平行，且小平台高度应与望远镜和平行光管下部齐平（为使得之后三棱镜放 置位置合适，能够折射平行光管发出的光） 粗调好后
○1 望远镜的调节：
望远镜内部构造如图 2
1. 接通图 1-10 下方小灯，从望远镜目镜观察图 2-6 处出现绿光，调节图 2-7 装
角
【实验数据与数据处理】
1. 实验条件：姓名：曾毅强 座位号：15 三棱镜编号：15 仪 =1‘ 入射光方位：10 = 3045’ ；20 = 1244’
2.实验数据 用自准法测量三棱镜定角 A 的数据记录表：
测量序号
第一位置 第二位置
1 游标Ⅰ 游标Ⅱ
740’ 18740’ 24745’ 6745’
○2 细调望远镜⊥分光计主轴 将小平台旋转 180°后捕捉到另一个 2-8 所示的像，但是此时该像并不在
清华大学分光计实验报告
大学物理实验报告分光计的调节和色散曲线的测定实验名称：分光计的调节和色散曲线的测定【实验目的】（1）了解分光计的原理与构造，学会调节分光计；（2）用最小偏向角法测定玻璃折射率；（3）掌握三棱镜顶角的两种测量方法。
一、分光计的调节及其原理（1）分光计的主要装置：望远镜，小平台，度盘，平行光管氦光谱管主要配件：三棱镜，平面镜分光计外形如图1（2）分光计的调节先进行粗调，通过自己的感觉，调整使得望远镜，小平台，平行光管与度盘平面平行，且小平台高度应与望远镜和平行光管下部齐平（为使得之后三棱镜放置位置合适，能够折射平行光管发出的光）粗调好后○1望远镜的调节：望远镜内部构造如图21.接通图1-10下方小灯，从望远镜目镜观察图2-6处出现绿光，调节图2-7装置使得能看清图2-3中“”形2.往小平台上放置一个平面镜正对望远镜物镜且与载物台调平螺钉中的两个的连线平行，微微调节图1-14找到较强绿光闪现位置，调节1-12装置使得该绿光出现在视野中央，调节1-9装置使得其成像清晰，如图2-8所示，再调1-12将其置于2-8所处位置并调整1-9使其消除视差。
此时望远镜已经适于观察平行光。
○2细调望远镜⊥分光计主轴将小平台旋转180°后捕捉到另一个2-8所示的像，但是此时该像并不在2-8所处位置（若还在该位置，则此时小平台和望远镜⊥分光计主轴），此时先调节小平台连线与平面镜平行的螺钉，使像与目标位置距离缩短与目标位置距离的一半，再通过调节1-9使其达到目标位置；再将小平台旋转180°，重复上述操作，直至像始终在2-8位置。
此时望远镜⊥分光计主轴。
○3调节平行光管1.调节平行光管发出平行光通过望远镜观察平行光管，调节1-28使得狭缝像应是一条轮廓清楚的窄长条，而不是边缘模糊的亮条；狭缝像应与叉丝无视差；狭缝宽度应合适，实验中缝像约１mm宽即可。
2.调节平行光管光轴⊥分光计主轴调节平行光管的调水平螺钉，使狭缝像被‘'叉丝中心线平分，平行光管也就垂直于分光计主轴了。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验六非线性方程求解实验目的1. 掌握用matlab软件求解非线性方程和方程组的基本用法, 并对结果做初步分析.2. 练习用非线性方程和方程组建立实际问题的模型并进行求解.实验内容题目3（1）小张夫妇以按揭方式贷款买了1 套价值20 万元的房子，首付了5 万元，每月还款1000 元，15 年还清。
问贷款利率是多少？（2）某人欲贷款50 万元购房，他咨询了两家银行，第一家银行开出的条件是每月还4500 元，15 年还清；第二家银行开出的条件是每年还450000 元，20 年还清。
从利率方面看，哪家银行较优惠（简单地假设年利率=月利率×12）？建立模型：设房价为b，首付款为b0，银行按照月利率(复利)来计算，月利率为r，月付款(月末支付)为a，共需要支付的月数为n。
根据经济学中资金的时间价值概念，可以得到：房价在n个月之后的实际价值为：b(1+r)n按揭购房期间交的所有款项在第n个月末的实际价值为：b0(1+r)n+a(1+r)n−1+(1+r)n−2+⋯+1=b0(1+r)n+a×(1+r)n−1由于在第n个月末还清了贷款，因此上述两个时间价值相等，则得到下面的关系式，即为解答此问题的方程：b(1+r)n=b0(1+r)n+a×(1+r)n−1即：(b−b0)(1+r)n−a×(1+r)n−1=0(1)代入已知条件：b=200000，b0=50000，a=1000，n=180，利用MATLAB解此非线性方程，经过简单的估测之后，给定初始值为r0=0.001，得到结果为：r=0.0020812，即贷款月利率为0.20812%。
（2）I.第一家银行相应的已知条件为：b=500000，b0=0，a=4500，n=180，利用MATLAB计算，经过简单的估测之后，给定初始值为r0=0.005，得到结果为：r=0.0058508，即这家银行的贷款月利率为0.58508%。
II.第二家银行由于按照年利率计算，因此方程中相应参数的意义有所改变，故已知条件为：b=500000，b0=0，a=45000，n=20，利用MATLAB计算，经过简单的估测之后，给定初始值也为r0=0.06，得到结果为。
r=0.063949，即这家银行的贷款年利率为6.3949%，则月利率为0.53291%实验结果：以月利息为比较条件，第二家银行比较优惠。
结果分析：(1)本题第二问里，将第二家银行的年利率近似看作是月利率的12倍，会造成一定的误差，若将计算结果0.53291%代入方程，可以得出每月需交付3697元，这样一年需要交付44364元，比题目中的45000元略小一些，说明计算得到的月利率偏小，但不至于影响最终的判断。
问题6：给定4种物质对应的参数ai, bi, ci和交互作用矩阵Q如下:a1=18.607, a2=15.841, a3=20.443, a4=19.293;b1=2643.31, b2=2755.64, b3=4628.96, b4=4117.07;c1=239.73, c2=219.16, c3=252.64, c4=227.44;Q=[ 1.0 0.192 2.169 1.6110.316 1.0 0.477 0.5240.377 0.360 1.0 0.2960.524 0.282 2.065 1.0]在压强p=760mmHg下，为了形成均相共沸混合物，温度和组分分别是多少？请尽量找出所有的可能解。
解: 设该混合物由n个可能的组分组成，组分i所占的比例为xi(i=1, … , n)，则∑_(i=1)^n▒〖x(i)〗=1, xi>=0 -------(1)xi((b(i))/(T+c(i))+ln(∑_(j=1)^n▒〖x(j)〗q(ij))+∑_(j=1)^n▒(x(j)q(ij))/(∑_(k=1)^n▒〖x(k)q(jk)〗) -1-aij+lnP)=0, i=1, … , n.-------(2)qij表示组分i与组分j的交互作用参数，qij构成交互作用矩阵Q程序:function f =azeofun(XT,n,P,a,b,c,Q)x(n)=1;for i=1:n-1x(i)=XT(i);x(n)=x(n)-x(i);endT=XT(n);p=log(P);for i=1:nd(i)=x*Q(i,1:n)';dd(i)=x(i)/d(i);endfor i=1:nf(i)=x(i)*(b(i)/(T+c(i))+log(x*Q(i,1:n)')+dd*Q(1:n,i)-a(i)-1+p);endn=4;P=760;a=[18.607, 15.841, 20.443, 19.293]';b=[2643.31, 2755.64, 4628.96, 4117.07];c=[239.73, 219.16, 252.64, 227.44];Q=[1.0 0.192 2.169 1.6110.316 1.0 0.477 0.5240.377 0.360 1.0 0.2960.524 0.282 2.065 1.0];XT0=[0.25,0.25,0.25,50];[XT,Y]=fsolve(@azeofun,XT0,[],n,P,a,b,c,Q)结果:XT=[0.0000 0.5858 0.4142 71.9657]Y= 1.0e-006 *[-0.0009 -0.0422 0.4428 -0.4701]分析:在上面计算中，对初值XT0的取法是:4种物质各占1/4,温度为50。
C。
初值解XT0 x1 x2 x3 x4 T[0.25,0.25,0.25,50] 0.0000 0.5858 0.4142 0.0000 71.9657[0.7,0.9,0.5,100] 0.0000 1.1425 0.0182 -0.1607 87.2356[0,1,0,72] 0.0000 0.7803 0.0000 0.2197 76.9613[0,0,1,72] -0.0000 0.0000 1.0000 0.0000 82.5567[0,0,0,90] -0.0000 -0.0000 -0.0000 1.0000 97.7712[1,0,0,30] 1.0000 -0.0000 -0.0000 0.0000 -18.9700………………………………分析:(i) 第一、三行解分别代表了组分二、三与组分二、四形成的均相共沸物，其组成与温度均符合要求；(ii)第四、五、六行解中均只有一种组分，违背了均相共沸物是由两种或两种以上物质组成的液体混合物的定义，故这些解均不符合要求；(iii)第二行第四种组分为负数，不符合要求。
即符合要求的均相共沸物的组成与温度分别是(i)第一种组分占0，第二种组分占0.5858，第三种组分占0.4142，最后一种组分占0，温度为71.9657℃；(ii)第一种组分占0，第二种组分占0.7803，第三种组分占0，最后一种组分占0.2197，温度为76、问题8：假设商品在t时刻的市场价格为p(t)，需求函数为D(p(t))=c-dp(t);(c,d>0)，而生产方的期望价格为q(t)，供应函数S(q(t))，当功效平衡是S(q(t))=D(p(t))。
若期望截个与市场价格不符，商品市场部均衡，生产方t+1时期的期望价格将会调整方式为q(t+1)-q(t)=r[p(t)-q(t)];(0<r<1)，以P(t)=[c-D(p(t))]/d=[c-S(q(t))]/d代入，得到关于q(t)的递推方程，设S(x)=arctan(ux)，u=4.8,d=0.25,r=0.3，以c为可变参数，讨论期望价格q(t)的变化规律，是否有混沌现象出现，并找出前几个分叉点，观察分叉点的极限趋势是否符合Feigenbaum常数揭示的规律。
模型建立：q的递推关系式如下：q(t+1)=q(t)+r[p(t)-q(t)]=q(t)+r*[c-arctan(uq(t))/d-q(t)]=(1-r)q(t)+r*c/d-r* arctan(uq(t))/d模型求解：建立chaos.m的源文件：function chaos(iter_fun,x0,r,n)kr=0;for rr=r(1):r(3):r(2)kr=kr+1;y(kr,1)=feval(iter_fun,x0,rr);for i=2:n(2)y(kr,i)=feval(iter_fun,y(kr,i-1),rr);endendplot([r(1):r(3):r(2)],y(:,n(1)+1:n(2)),'k.');建立iter01.m的源文件：function y=iter01(x,c)u=4.8,d=0.25,r=0.3;y=(1-r)*x+r*c/d-r*atan(u*x)/d;主程序：chaos(@iter01,0.5,[0.2,1.6,0.01],[100,200])得到的结果如下：改变chaos.m中的坐标，我们可得到相应分叉点为：1.086 0.953 0.907 0.897 (b2-b1)/(b3-b2)=2.8913(b2-b2)/(b4-b3)=4.6000由此得出，当n越大时比值越接近于Feigenbaum常数：4.6692，因此满足此规律。
模型检验：建立iter02.m的源文件：function y=iter02(c)u=4.8;d=0.25;r=0.3;q(1)=0.5;for n=1:1:50q(n+1)=(1-r)*q(n)-r*atan(u*q(n))/d+r*c/d;endN=1:1:51;plot(N,q);主程序如下：>> figure;iter02(1.1);>> figure;iter02(1);>> figure;iter02(0.93);>> figure;iter02(0.9);>> figure;iter02(0.7);得到的结果依次如下：C=1.1，模型收敛于一点当c=1：模型在两点之间震荡。
当c=0.93：模型在4点之间震荡。
当c=0.9：模型在8点之间震荡，呈现初始混沌状态。
当c=0.7：模型处于混沌状态。
从上述实验结果得出混沌现象的出现与分叉点十分符合，分叉点是极限趋势是符合Feigenbaum常数揭示的规律。