清华大学数学实验报告4
- 格式:docx
- 大小:207.83 KB
- 文档页数:25
第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. 理解线性方程组的概念和求解方法。
2. 掌握高斯消元法和矩阵求逆法求解线性方程组。
3. 熟悉MATLAB软件在数学实验中的应用。
实验时间:2021年X月X日实验地点:计算机实验室实验器材:1. 计算机2. MATLAB软件实验内容:一、实验原理线性方程组是数学中一类常见的方程组,其形式如下:\[ Ax = b \]其中,\( A \) 是一个 \( m \times n \) 的系数矩阵,\( x \) 是一个 \( n \) 维的未知向量,\( b \) 是一个 \( m \) 维的常数向量。
线性方程组的求解方法有多种,如高斯消元法、矩阵求逆法等。
本实验主要介绍高斯消元法和矩阵求逆法。
二、实验步骤1. 设计一个线性方程组,并记录系数矩阵 \( A \) 和常数向量 \( b \)。
\[ \begin{cases}2x + 3y - z = 8 \\-x + 2y + 3z = 1 \\4x - y + 2z = 3\end{cases} \]系数矩阵 \( A \) 和常数向量 \( b \) 如下:\[ A = \begin{bmatrix}2 &3 & -1 \\-1 & 2 & 3 \\4 & -1 & 2\end{bmatrix}, \quad b = \begin{bmatrix}8 \\1 \\3\end{bmatrix} \]2. 使用MATLAB软件进行高斯消元法求解线性方程组。
```matlabA = [2 3 -1; -1 2 3; 4 -1 2];b = [8; 1; 3];x = A\b;```3. 使用MATLAB软件进行矩阵求逆法求解线性方程组。
```matlabA_inv = inv(A);x_inv = A_invb;```4. 比较两种方法得到的解,并验证其正确性。
三、实验结果与分析1. 使用高斯消元法求解得到的解为:\[ x = \begin{bmatrix}2 \\1 \\1\end{bmatrix} \]2. 使用矩阵求逆法求解得到的解为:\[ x = \begin{bmatrix}2 \\1 \\1\end{bmatrix} \]两种方法得到的解相同,验证了实验的正确性。
回归分析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做二元线性回归。
数学实验报告样本标题:投影性质实验报告一、引言投影是数学中一个重要的概念,它在几何学、线性代数以及物理学等领域中都有广泛的应用。
本实验旨在通过实际操作和观察,探究几何图形在不同投影方式下的性质。
二、实验内容1.准备材料:白色纸张、直尺、铅笔、胶带。
2.实验步骤:a.在纸张上画出一些几何图形,如三角形、矩形、正方形等。
b.选择一个固定点作为观察点,将纸张用胶带固定在观察点上方。
c.将光源放置在观察点的正后方,以确保光线垂直投射到纸张上。
d.观察并记录图形在纸张上的投影。
三、实验结果1.绘制图形:我们选择绘制了一个三角形、一个矩形和一个正方形作为实验对象,并将它们固定在观察点上方。
这样可以保证光线从正上方垂直投射到纸上的每个图形。
2.观察结果:a.三角形的投影是一个三角形,其形状与原图形相似,但是大小可能会有所不同。
b.矩形的投影是一个矩形,其形状与原图形相同。
c.正方形的投影是一个正方形,其形状与原图形相同。
3.结果分析:从观察结果可以看出,当几何图形与观察点和光源的位置关系较为简单时,其投影形状与原图形相似。
特别是在观察点和光源位置固定的情况下,图形的大小可能会有所改变,但形状保持不变。
四、讨论1.关于投影形状:每种几何图形在不同的投影方式下可能会有不同的形状。
投影形状的变化取决于观察点和光源的位置关系、以及几何图形本身的性质。
2.关于投影大小:在本实验中,我们观察到图形的大小可能会发生变化。
这是由于观察点和光源的位置决定了图形在纸上的投影长度。
当观察点与光源距离增加时,投影相对于原图形可能会变大;反之,当距离减少时,投影可能会变小。
3.关于应用:投影性质是计算机图形学、建筑设计以及摄影学等领域中的关键概念之一、准确理解和运用投影性质可以帮助我们更好地设计和呈现图形。
五、结论通过本实验,我们实际操作和观察了几何图形在不同投影方式下的性质。
我们观察到,在固定观察点和光源位置的情况下,图形的形状保持不变,但大小可能会发生变化。
一、实验目的1. 理解信号的基本运算概念,包括信号的加法、减法、乘法和除法。
2. 掌握使用MATLAB进行信号运算的方法。
3. 分析信号运算后的特性,如幅度、相位和时域变化。
二、实验原理信号的运算是指对两个或多个信号进行数学运算,得到新的信号。
常见的信号运算包括:1. 信号的加法:将两个信号的幅度值相加,得到新的信号。
2. 信号的减法:将一个信号的幅度值减去另一个信号的幅度值,得到新的信号。
3. 信号的乘法:将两个信号的幅度值相乘,得到新的信号。
4. 信号的除法:将一个信号的幅度值除以另一个信号的幅度值,得到新的信号。
三、实验仪器与软件1. 仪器:示波器、信号发生器、计算机2. 软件:MATLAB四、实验内容与步骤1. 实验一:信号的加法与减法(1)使用信号发生器产生两个正弦信号,频率分别为1Hz和2Hz,幅度分别为1V和2V。
(2)将两个信号分别输入示波器,观察波形。
(3)使用MATLAB编写程序,将两个信号相加和相减,并绘制结果波形。
(4)分析结果,比较加法和减法运算对信号特性的影响。
2. 实验二:信号的乘法与除法(1)使用信号发生器产生两个正弦信号,频率分别为1Hz和2Hz,幅度分别为1V和2V。
(2)将两个信号分别输入示波器,观察波形。
(3)使用MATLAB编写程序,将两个信号相乘和相除,并绘制结果波形。
(4)分析结果,比较乘法和除法运算对信号特性的影响。
3. 实验三:信号运算的时域分析(1)使用MATLAB编写程序,对实验一和实验二中的信号进行时域分析,包括信号的幅度、相位和时域变化。
(2)比较不同信号运算后的特性变化。
五、实验结果与分析1. 实验一:信号的加法与减法通过实验,观察到信号的加法和减法运算对信号的幅度和相位有显著影响。
加法运算使信号的幅度增加,相位保持不变;减法运算使信号的幅度减小,相位保持不变。
2. 实验二:信号的乘法与除法通过实验,观察到信号的乘法和除法运算对信号的幅度和相位有显著影响。
一、实验目的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\) 随时间的变化呈现出无规则、复杂的振荡行为。
实验4 液体饱和蒸气压的测定化53卢巍2005011871 实验同组人:须丹丹 实验日期:071030 交实验报告日期:071113 实验指导教师姓名:贾维杰1、 引言:1.1实验目的:1. 运用克劳修斯-克拉贝龙方程,求出所测温度范围内的平均摩尔气化焓及正常沸点。
2. 掌握测定饱和蒸汽压的方法。
1.2实验原理在通常温度下(距离临界温度较远时),纯液体与其蒸气达平衡时的蒸气压称为该温度下液体的饱和蒸气压,简称为蒸气压。
蒸发1摩尔液体所吸收的热量称为该温度下液体的摩尔气化热。
液体的蒸气压与液体的本性及温度等因素有关。
随温度不同而变化,温度升高时,蒸气压增大;温度降低时,蒸气压降低,这主要与分子的动能有关。
当蒸气压等于外界压力时,液体便沸腾,此时的温度称为沸点,外压不同时,液体沸点将相应改变,当外压为p ø(101.325kPa )时,液体的沸点称为该液体的正常沸点。
液体的饱和蒸气压与温度的关系用克劳修斯(Clausius )-克拉贝龙(Clapeyron )方程式表示:式中,R 为摩尔气体常数;T 为热力学温度;Δvap H m 为在温度T 时纯液体的摩尔气化热。
假定Δvap H m 与温度无关,或因温度变化范围较小,Δvap H m 可以近似作为常数,积分上式,得:Aln p B T=-+ 或 Aln p B T=-+ vap m H Rm ∆=-式中:B ——积分常数。
从上式可知:若将ln p 对1/T 作图应得一直线,斜率m=vap m -A H /R =-∆。
由此可得 vap m H Rm ∆=-同时从图上可求出标准压力时的正常沸点。
2、实验操作:2.1 实验药品、仪器型号及测试装置示意图实验药品及仪器:等压管1支、稳压瓶1个、负压瓶1个、恒温槽1套、真空泵1台、压力计1台。
乙醇(分析纯)图2-4-1 纯液体饱和蒸气压测定装置:1、等压管,2、冷凝管,3、搅拌器,4、加热器,5、1/10゜C温度计,6、辅助温度计,7、稳压瓶,8、负压瓶,9、干燥管。
数学实验内容及目的(作业参考)实验1 一元函数的图形实验目的1.学习 matlab一元函数绘图命令.2.进一步理解函数概念.实验内容1.学习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 可打开“函数计算器”图形用户界面.2.绘制函数图形举例. 例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 函数极限实验目的1.理解极限概念.2.掌握用matlab 软件求函数极限的方法.实验内容1.学习matlab 命令.matlab 求极限命令可列表如下:表2.1matlab 代数方程求解命令solve 调用格式. solve (函数)(x f ) 给出0)(=x f 的根.2.理解极限概念.数列}{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 =1.对函数的极限概念,也可用上述方法理解.例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之间无限次振荡,极限不存在.仔细观察该图象,发现图象的某些峰值不是1和-1,而我们知道正弦曲线的峰值是1和-1,这是由于自变量的数据点选取未必使x 1sin 取到1和-1的缘故,读者可试增加数据点,比较它们的结果. 例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 . 3.求函数极限例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)4.求方程的解. 例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 导数及偏导数计算实验目的1.进一步理解导数概念及其几何意义. 2.学习matlab 的求导命令与求导法.实验内容1.学习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 f2.导数概念.导数是函数的变化率,几何意义是曲线在一点处的切线斜率. (1)点导数是一个极限值.例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.可知1)0(='f(2)导数的几何意义是曲线的切线斜率.例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 越来越接近曲线的割线. 3.求一元函数的导数.(1))(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 个分量.(2)参数方程所确定的函数的导数.设参数方程⎩⎨⎧==)()(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)). 其中分号的作用是不显示结果. 4.求多元函数的偏导数.例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)]. 5.求高阶导数或高阶偏导数.例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.6.求隐函数所确定函数的导数或偏导数 例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 积分计算实验目的1.通过本实验加深理解积分理论中分割、近似、求和、取极限的思想方法. 2.学习并掌握用matlab 求不定积分、定积分、二重积分、曲线积分的方法. 3.学习matlab 命令sum 、symsum 与int . 实验内容1.学习matlab 命令(1)求和命令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(2)求和命令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 (3)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 ),(2.计算不定积分 例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)] 3.定积分的概念.定积分是一个和的极限.取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)-14.计算定积分和广义积分.例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)的值.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.6.曲线积分例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 自定义函数与导数应用实验目的1.学习matlab 自定义函数.2.加深理解罗必塔法则、极值、最值、单调性. 实验内容1.学习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 上的最小值点. 2.自定义函数例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.76562.验证罗必塔法则.罗必塔法则是指在求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)(lim4.函数的单调性与极值.例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 . 5.函数的最值调用求函数最小值命令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 矩阵运算与数组运算实验目的:1.理解矩阵及数组概念.2.掌握matlab 对矩阵及数组的操作命令.实验内容:1.矩阵与数组的输入. 对于较小较简单的矩阵,从键盘上直接输入矩阵是最常用的数值矩阵创建方法.用这种方法输入矩阵时注意以下三点:(1)整个输入矩阵以方括号“[ ]”为其首尾; (2)矩阵的元素必须以逗号“,”或空格分隔; (3)矩阵的行与行之间必须用分号“;”或回车键隔离. 例1:下面的指令可以建立一个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列的矩阵,数组的输入可以采用上面矩阵的输入方法.例2:输入以下指令以建立数组c.c=[1 2 3 4 5 6 7 8]c =1 2 3 4 5 6 7 8另外还有两种方法输入数组.请看下面两个例子.例3:在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注意“:”的使用方法.例4:在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个数据点的数组.2.常用矩阵的生成.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的行数和列数.例5:生成以下矩阵.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 03.矩阵元素的标识矩阵的元素、子矩阵可以通过标量、向量、冒号的标识来援引和赋值.(1)矩阵元素的标识方式A(ni,nj).ni,nj都是标量.若它们不是整数,则在调用格式中会自动圆整到最临近整数.ni指定元素的行位置,nj指定元素的列位置.(2)子矩阵的序号向量标识方式A(v,w).v,w是向量,v,w中的任意一个可以是冒号“:”,他表示取全部行(在v位置)或全部列(在w位置).v,w中所用序号必须大于等于1且小于等于矩阵的行列数.例6:元素和矩阵的标识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 04.矩阵运算和数组运算.矩阵运算的指令和意义如下: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相乘.例7: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的元素为指数求幂.例8: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矩阵与线性方程组实验目的:1.掌握matlab求矩阵的秩命令.2.掌握matlab求方阵的行列式命令.3.理解逆矩阵概念,掌握matlab求逆矩阵命令.4.会用matlab求解线性方程组.实验内容:1.矩阵的秩.指令rank(A)将给出矩阵A的秩.例1: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 =22.方阵的行列式.指令det(A)给出方阵A的行列式.例2: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你能解释上例中的运算结果吗?在这里我们实际上验证了行列式的性质.3.逆矩阵指令inv(A)给出方阵A的逆矩阵,如果A不可逆,则inv(A)给出的矩阵的元素都是Inf.例3:设⎪⎪⎪⎭⎫⎝⎛=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 <时,矩阵除给出方程的最小范数解.例4:解方程组:⎪⎪⎩⎪⎪⎨⎧=-+=++=+-+=++-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例5:解方程组:⎪⎩⎪⎨⎧=-++-=-++-=--++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 +=例6:解方程组:⎪⎪⎩⎪⎪⎨⎧=+--=--=-+-=+--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 常微分方程与级数实验目的:1.学习用matlab求解微分方程命令dsolve.2.学习matlab泰勒级数展开命令.3.巩固幂级数的收敛半径、和等概念.实验内容:1.学习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,自变量的趋向值)2.求解一阶微分方程.微分方程在输入时,y'应输入Dy,y''应输入D2y等,D应大写.例1:求微分方程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)例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))例3:求微分方程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)3.求解二阶微分方程.例4:求03=+'+''x e y y 的通解.解:输入命令:dsolve('D2y+3*Dy+exp(x)=0','x')得结果:ans =-1/4*exp(x)+C1+C2*exp(-3*x)例5:求解微分方程.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*C14.taylor 展开式.例6:求函数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 5.级数求和.例7:求∑∞=121n n .解:输入命令:syms n;symsum(1/2^n,1,inf) 得结果:ans =1例8:求幂级数∑∞=⨯12nnnnx的和函数.解:输入命令:symsum(x^n/(n*2^n),n,1,inf) 得结果ans =-log(1-1/2*x)例9:求幂级数∑∞=1nnnx的和函数.解:输入命令:symsum(n*x^n,n,1,inf) 得结果:ans =x/(x-1)^26.判别级数敛散性.例10:判断数项级数∑∞=+1)1(1nnn的收敛性.解:输入求和命令:symsum(1/(n*(n+1)),n,1,inf) 得结果:ans =1求和得是1,说明该级数收敛.例11:判别级数∑∞=+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级数收敛,得知所要判别的级数也收敛.例12:判别级数∑∞=⎪⎭⎫⎝⎛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。
实验4 MATLAB 程序设计一、实验目的:1.理解命令式M 文件和函数式M 文件的关系及各自的特点;2.掌握MATLAB 中关系运算和逻辑运算的正确使用方法;3.掌握MATLAB 中两种循环结构及选择结构的正确使用方法;4.能够编制简单的程序解决实际问题。
二、实验内容:1.满足条件01=a ,12=a ,…,n n n a a a +=++12(n=0,1,…)的数列称为Fibonacci 数列,(1)编写一个MATLAB 函数式M 文件fibo.m 来计算Fibonacci 数列的任意项的数值;(2)计算数列中第50项的数值并给出程序运行所用的时间;(3)编写一个MATLAB 命令式M 文件调用函数式M 文件fibo.m 来计算Fibonacci 数列的前n 项和的值;给出n=50时的结果。
(1)计算Fibonacci 数列的程序function x=fibo(n)if n<1|fix(n)~=ndisp('N must be a positive integer.');elseif n==1;x=0;elseif n==2x=[0,1];elsex=zeros(1,n);x(1)=0;x(2)=1;for k=3:nx(k)=x(k-1)+x(k-2);endend(2)计算数列中第50项的数值x= fibo(50);ticF=x(50)toc(3)Fibonacci 数列的前50项和的值x= fibo(50);FF=sum(x)2.利用下面的公式计算π的值14134171513114---++-+-=N N π,选择不同的N 值,比较结果的精确度。
function PI=hw2(n)if n<1|fix(n)~=ndisp('N must be a positive integer.');elsePI=0;for k=1:nx=1/(4*k-3)-1/(4*k-1);PI=PI+x;endendPI=4*PI;n N的最小正整数n。
清华大学数学实验报告4————————————————————————————————作者: ————————————————————————————————日期:ﻩ电13 苗键强2011010645一、实验目的1.掌握用 MATLAB 软件求解非线性方程和方程组的基本用法, 并对结果作初步分析;2.练习用非线性方程和方程组建立实际问题的模型并进行求解。
二、实验内容题目1【问题描述】(Q1)小张夫妇以按揭方式贷款买了1套价值20万元的房子,首付了5万元,每月还款1000元,15年还清。
问贷款利率是多少?(Q2)某人欲贷款50 万元购房,他咨询了两家银行,第一家银行开出的条件是每月还4500元,15 年还清;第二家银行开出的条件是每年还45000 元,20年还清。
从利率方面看,哪家银行较优惠(简单假设:年利率=月利率×12)?【分析与解】假设初始贷款金额为x0,贷款利率为p,每月还款金额为x,第i 个月还完当月贷款后所欠银行的金额为x i,(i=1,2,3,......,n)。
由题意可知:x1=x0(1+p)−xx2=x0(1+p)2−x(1+p)−xx3=x0(1+p)3−x(1+p)2−x(1+p)−x……x n=x0(1+p)n−x(1+p)n−1−⋯−x(1+p)−x=x0(1+p)n−x (1+p)n−1p=0因而有:x0(1+p)n=x (1+p)n−1p (1)则可以根据上述方程描述的函数关系求解相应的变量。
(Q1)根据公式(1),可以得到以下方程:150p(1+p)180−(1+p)180+1=0设 f(p)=150p(1+p)180−(1+p)180+1,通过计算机程序绘制f(p)的图像以判断解p的大致区间,在Matlab中编程如下:fori = 1:25t = 0.0001*i;p(i) = t;f(i) =150*t*(1+t).^180-(1+t).^180+1;end;plot(p,f),hold on,grid on;运行以上代码得到如下图像:f(p)~p关系曲线图通过观察上图可知p∈[0.002,0.0022]。
Solution1:对于p∈[0.002,0.0022],采用二分法求解,在Matla b中编程如下:clear;clc;x0=150000;n =180;x = 1000;p0 = 0.002;p1=0.0022;while (abs(p1-p0)>1e-8)f0 = x0*(1+p0).^n+x*(1-(1+p0).^n)/p0;f1= x0*(1+p1).^n+x*(1-(1+p1).^n)/p1;p2=(p0+p1)/2;f2 = x0*(1+p2).^n+x*(1-(1+p2).^n)/p2;if (f0*f2>0 && f1*f2<0)p0= p2;elsep1 = p2;end;end;p0结果得到p0=0.078125=0.2081%.所以贷款利率是0.2081%。
Solution2:对于p∈[0.002,0.0022],采用牛顿法求解,为观测p{k}是否收敛,在Matlab中编程如下:clearclcn = 5;fori= 1:np(i) = 0.0001*(i+18);t =p(i);f(i) = t-(150*t*(1+t).^180-(1+t).^180+1)/(27000*t*(1+t).^179+150*(1+t).^180-180*(1+t).^179);g(i) = t;end;plot(p,f,p,g),hold on,grid on;运行以上代码得到如下图像:收由图像可知蓝色曲线在两线交点处斜率绝对值小于1,故p{k}敛。
取初始值p=0.0019,采用牛顿法求解,在Matlab中编程如下: p= 0.0019;for i=1:100a = p-(150*p*(1+p).^180-(1+p).^180+1)/(27000*p*(1+p).^179+150*(1+p).^180-180*(1+p).^179);p=a;endp结果得到p0=0.945915=0.2081%.所以贷款利率是0.2081%。
Solution3:采用fzero求解,在Matlab中编程如下:p= fzero(inline('150*x*(1+x).^180+1-(1+x).^180'),[0.00200,0.00225]);结果得到p0=0.945920=0.2081%.所以贷款利率是0.2081%。
【结论】贷款利率是0.2081%。
(Q2)根据公式(1),对于第一家银行提供的条件可以得到关于月利率p1的方程:1000p1(1+p1)180−9(1+p1)180+9=0对于第二家银行提供的条件可以得到关于年利率p2的方程:100p2(1+p2)20−9(1+p2)20+9=0设 f(p1)=1000p1(1+p1)180−9(1+p1)180+9,f(p2)=100p2(1+p2)20−9(1+p2)20+9。
通过计算机程序绘制f(p1)和f(p2)的图像以判断解p1和p2的大致区间,在Matlab中编程如下:clearclcfor i =1:100p(i) = 0.0001*i;t = p(i);f1(i) =1000*t*(1+t).^180-9*(1+t).^180+9;end;plot(p,f1),hold on,grid on;for i = 1:1000p(i) = 0.0001*i;t= p(i);f2(i) = 100*t*(1+t).^20-9*(1+t).^20+9;end;plot(p,f2),hold on,grid on;运行以上代码得到如下图像:f(p1)~p1关系曲线图f(p2)~p2关系曲线图通过观察以上图像可知p1∈[0.0055,0.0060], p2∈[0.060,0.065] Solution1:采用二分法求解相应的利率值,在Matlab中编程如下:clearclcp10 =0.0055;p11 = 0.0060;while(abs(p10-p11)>1e-8)f0 = 1000*p10*(1+p10).^180-9*(1+p10).^180+9;f1 = 1000*p11*(1+p11).^180-9*(1+p11).^180+9;p22=(p10+p11)/2;f2 =1000*p22*(1+p22).^180-9*(1+p22).^180+9;if (f0*f2>0 && f1*f2<0)p10 = p22;elsep11 = p22;end;end;p10clearclcp20=0.060;p21= 0.065;while (abs(p20-p21)>1e-8)f0= 100*p20*(1+p20).^20-9*(1+p20).^20+9;f1 = 100*p21*(1+p21).^20-9*(1+p21).^20+9;p22 = (p20+p21)/2;f2=100*p22*(1+p22).^20-9*(1+p22).^20+9;if (f0*f2>0 && f1*f2<0)p20 = p22;elsep21 = p22;end;end;p20结果得到p1=0.115234=0.5851%,相应的年利率为p1n=12*p1=7.021%;p2=0.77685=6.395%。
所以第一家银行的贷款年利率是7.021%,第二家银行的贷款年利率是6.395%。
Solution2:采用牛顿法求解相应的利率值,为观察p1{k}和p2{k}是否收敛,在Matlab中编程如下:clearclcn=6;fori=1:np(i)=0.0001*(i+55);t=p(i);f1(i)=t-(1000*t*(1+t).^180-9*(1+t).^180+9)/(180000*t*(1+t).^179+1000*(1+t).^180-1620*(1+t).^179);g(i)=t;end;plot(p,f1,p,g),hold on,grid on;clearclcn = 10;fori= 1:np(i) = 0.001*(i+60);t= p(i);f2(i)= t-(100*t*(1+t).^20-9* (1+t).^20+9)/(2000*t*(1+t).^19+100*(1+t).^20-180*(1+t).^19);g(i)=t;end;plot(p,f2,p,g),holdon,grid on;运行以上代码得到如下图像:由图像可知两条蓝色曲线在两线交点处斜率绝对值均小于1,故p1{k}和p2{k}均收敛。
,取初始值p=0.0056,对于p2{k},取初始值p=0.061,对于p1{k}采用牛顿法求解,在Matlab中编程如下:clearclcp =0.0056;for i = 1:100t =p;pp = t-(1000*t*(1+t).^180-9*(1+t).^180+9)/(180000*t *(1+t).^179 +1000*(1+t).^180-1620*(1+t).^179);p= pp;endpclearclcp= 0.061;for i = 1:100t = p;pp= t-(100*t*(1+t).^20-9*(1+t).^20+9)/(2000*t*(1+t).^19+100*(1+t).^20-180*(1+t).^19);p= pp;endp结果得到p1=0.284539=0.5851%,相应的年利率为p1n=12*p1=7.021%;p2=0.23863=6.395%。
所以第一家银行的贷款年利率是7.021%,第二家银行的贷款年利率是6.395%。
Solution3:采用fzero求解,在Matlab中编程如下:clearclcp1 = fzero(inline('1000*x*(1+x).^180+9-9*(1+x).^180'),[0.0056,0.0066]);p2 = fzero(inline('100*x*(1+x).^20+9-9*(1+x).^20'), [0.061,0.071]);p1p2结果得到p1=0.284532=0.5851%,相应的年利率为p1n=12*p1=7.021%;p2=0.23863=6.395%。
所以第一家银行的贷款年利率是7.021%,第二家银行的贷款年利率是6.395%。
【结论】由于第一家银行的贷款年利率是7.021%,第二家银行的贷款年利率是6.395%,所以第二家银行较优惠。