用matlab解根3乘根2,MATLAB应用.ppt

资源描述

数学建模应用软件培训 MATLAB应用哈尔滨理工大学应用数学系 1 方程求根2 线性方程组求解3 插值 拟合4 应用 人口问题 划艇比赛成绩 主要内容 1 方程求根2 线性方程组求解3 插值 拟合4 应用 人口问题 划艇比赛成绩 主要内容 1方程求根 1 1方程求根1 2MATLAB主要命令介绍 1 1方程求根 问题 求方程f x 0的解方法 因式分解法 图形放大法 数值迭代逼近法 迭代逼近法 区间的迭代 点的迭代区间迭代又分为对分法和黄金分割法 点的迭代又分为简单迭代法 单点割线法 两点割线法 牛顿法等 1 2MATLAB主要命令介绍 多项式求根 求满足多项式p x 0的x值 N次多项式应该有n个根 这些根可能是实根 也可能是若干对共轭复根 调用格式 x roots P 其中P为多项式的系数向量 求得的根赋给向量x 即x 1 x 2 x n 分别代表多项式的n个根 注 该命令每次只能求一个一元多项式的根 该指令不能用于求方程组的解 必须把多项式方程变成Pn x 0的形式 1 2MATLAB主要命令介绍 例 求方程的解 首先将方程变成Pn x 0的形式 输入命令 roots 1 10 1 例 求多项式方程x4 8x3 10 0的根 A 1 8 0 0 10 x roots A x 8 01941 0344 0 5075 0 9736i 0 5075 0 9736i ans 1 4656 0 2328 0 7926i 0 2328 0 7926i 1 2MATLAB主要命令介绍 solve s 求解符号表达式s的代数方程 求解变量为默认变量 当方程右端为0时 方程可以不标出等号和0 仅标出方程的左端 solve s v 求解符号表达式s的代数方程 求解变量为v solve s1 s2 sn v1 v2 vn 求解符号表达式s1 s2 sn组成的代数方程组 求解变量分别v1 v2 vn solve 1 2MATLAB主要命令介绍 例1 解下列方程 1 x solve 1 x 2 4 x x 2 4 1 2 x 2 x 2 f sym x x 3 4 x 7 1 3 1 x solve f solve x 1 f x x 3 4 x 7 1 3 1x 3 1 2MATLAB主要命令介绍 3 x solve 2 sin 3 x pi 4 1 4 x solve x x exp x 10 x 仅标出方程的左端 solve x 5 pi 36 13 pi 36 x 1 6335061701558463841931651789789 1 2MATLAB主要命令介绍 求解方程f x 0的实数根也就是求函数f x 的零点 MATLAB中设有求函数f x 零点的指令fzero 可用它来求方程的实数根 fzero fun x0 options 输入参数fun为函数f x 的字符表达式 内联函数名或M函数文件名 输入参数x0为函数某个零点的大概位置 不要取零 或存在的区间 xi xj 要求函数f x 在x0点左右变号 即f xi f xj 0 函数零点指令fzero 1 2MATLAB主要命令介绍 输入参数options可有多种选择 若用optimset disp iter 代替options时 将输出寻找零点的中间数据 该指令无论对多项式函数还是超越函数都可以使用 但是每次只能求出函数的一个零点 因此在使用前需摸清函数零点数目和存在的大体范围 为此 一般先用绘图指令plot fplot或ezplot画出函数f x 的曲线 从图上估计出函数零点的位置 函数零点指令fzero 1 2MATLAB主要命令介绍 例 求方程x2 4sin x 25的实数根 2 x 2 解1 fun为函数f x 的字符表达式 l 首先要确定方程实数根存在的大致范围 为此 先将方程变成标准形式f x x2 4sin x 25 0 作f x 的曲线图 x 2 pi 0 1 2 pi f x 2 4 sin x 25 F是一个向量 点不能去掉 plot x f gridon 描绘点图 从曲线上可以看出 函数的零点大约在x1 4和x2 5附近 1 2MATLAB主要命令介绍 解2 fun为函数f x 的M函数文件名将方程x2 4sin x 25编成M函数文件 实用中的函数较为复杂 而又多次重复调用时 才这样做 用fzero求解 1 在M文件编辑调试窗中键入 functionyy li4 5 x yy x 2 4 sin x 25 以li4 5为文件名存盘 退出编辑调试窗 回到指令窗 1 2MATLAB主要命令介绍 解2 fun为函数f x 的M函数文件名 2 确定根的大体位置 3 在指令窗中键入下述指令可求出 4附近的根x1 fzero li4 5 4 键入下述指令可求出5附近的根 x2 fzero li4 5 5 1 2MATLAB主要命令介绍 2 直接使用指令fzero求出方程在x1 4时的根x1 fzero x 2 4 sin x 25 4 若键入 fzero x 2 4 sin x 25 4 optimset disp iter 将显示迭代过程 中间数据表明 求根过程中不断缩小探测范围 最后得出 4附近满足精度的近似根 3 求x2 5的根 x2 fzero x 2 4 sin x 25 5 1 2MATLAB主要命令介绍 符号方程组求解 在MATLAB中 求解用符号表达式表示的方程组仍然可由函数solve实现 其调用格式与解用符号表达式表示的方程一样 例 解 xy solve 1 x 3 1 y 3 28 1 x 1 y 4 x y x 11 3y 1 31 1 2MATLAB主要命令介绍 求方程组数值解的指令最小二乘解不是精确解 fsolve是用最小二乘法求解非线性方程组F X 0的指令 变量X可以是向量或矩阵 方程组可以由代数方程或者超越方程构成 它的使用格式为 fsolve fun X0 哪个点附近的零点 OPTIONS 参数fun是编辑并存盘的M函数文件的名称 可以用 代替单引号对它进行标识 M函数文件主要内容是方程F X 0中的函数F X 即方程左边的函数 1 2MATLAB主要命令介绍 参数X0是向量或矩阵 为探索方程组解的起始点 求解将从X0出发 逐渐趋向 最终得到满足精度要求 最接近X0的近似根X F X 0 由于X0是向量或矩阵 无法用画图方法进行估计 实际问题中常常是根据专业知识 物理意义等进行估计 该指令输出一个与X0同维的向量或矩阵 为方程组的近似数值解 参数OPTIONS为设置选项 用它可以设置过程显示与否 误差 算法 具体内容可用help查阅 通常可以省略该项内容 1 2MATLAB主要命令介绍 例求方程组在x0 1 y0 1 z0 1附近的数值解 解 l 在文本编辑调试窗中编辑M函数文件 首先将方程组变换成F X 0的形式 x y z看成向量X的三个分量 1 2MATLAB主要命令介绍 functionms li4 7 X ms 1 X 1 2 10 X 1 X 2 2 X 3 7 ms 2 X 1 X 2 2 2 X 3 ms 3 X 1 2 X 2 2 3 X 2 X 3 2 这里 X l x X 2 y X 3 z 输出参量ms也有三个分量用 li4 7 为M函数文件名存盘 1 2MATLAB主要命令介绍 2 在指令窗中键入 fsolve li4 7 111 若键入 x fsolve li4 7 111 optimset Display iter 则得出求解过程 该方程也可以用MATLAB的符号指令solve求解 但结果非常冗长 1 方程求根2 线性方程组求解3 插值 拟合4 应用 人口问题 划艇比赛成绩 主要内容 2线性方程组求解 2 1MATLAB主要命令介绍2 2线性方程组求解问题2 3应用实例 2 1MATLAB主要命令介绍 A B命令 A BA B A BAX B 当A可逆时 X A B 此时用左除 XA B 当A可逆时 X B A 此时用右除 1 矩阵加法 2 矩阵乘积 3 矩阵左除和右除 X A B X B Adet A det A A A inv A inv A 2 1MATLAB主要命令介绍 4 求矩阵的行列式 5 矩阵的转置 6 矩阵的逆矩阵 rank A rank A rref A rref A formatratdiag A 2 1MATLAB主要命令介绍 7 矩阵的秩 8 将矩阵化为最简行阶梯形矩阵 9 输出格式 10 对角阵 例2 1已知 求A B A 23 57 B 37 1155 A B A B 1101622 2 1MATLAB主要命令介绍 例2 2已知 求A B 解 A 23 57 B 375 11557 A B 输出 A B 2759316214074 2 1MATLAB主要命令介绍 例2 3已知 求X 解 输入 A 21 51 1 30 6 02 12 14 76 B 8 9 5 0 X A B 左除X 3 0000 4 0000 1 00001 0000 2 1MATLAB主要命令介绍 例2 4已知 求A的行列式 解 A 12 46 det A ans 2 2 1MATLAB主要命令介绍 例2 5求的行列式 解 symsabcd 符号变量必须要定义 A ab cd det A det A a d b c 2 1MATLAB主要命令介绍 A的行列式 创建矩阵 定义变量abcd 例2 6已知 求A的转置 解 A 15 15 11 3 A A 1 11155 3 2 1MATLAB主要命令介绍 例2 7将化为行最简形阶梯矩阵 解 A 12 46 rref A 输出 rref A 1 00000 0 861101 0000 0 8889000 2 1MATLAB主要命令介绍 2 1MATLAB主要命令介绍 如果在输入rref A 前 输入 formatrat rref A 输出 rref A 10 31 301 8 9000恢复原来的小数输出 输入 format 恢复输出原来的格式 例2 8创建以1 2 3为主对角线元素的对角矩阵 解 V 123 s diag V 或 s diag V 0 s 100020003 2 1MATLAB主要命令介绍 方程组 2 1 记则方程组 2 1 可表示为 2 2线性方程组求解 记即方程组 2 1 的增广矩阵 2 2线性方程组求解 齐次线性方程组解的讨论 一定有解 至少零解 且当时 有唯一零解 当时 有非零解 且有个线性无关的解向量 2 2线性方程组求解 非齐次线性方程组解的讨论 如果系数矩阵的秩不等于增广矩阵的秩 则无解 如果系数矩阵的秩等于增广矩阵的秩时 当秩为时 有唯一解 如果系数矩阵的秩等于增广矩阵的秩时 当秩小于时 有无穷多解 可求出通解 2 2线性方程组求解 超定方程设线性方程组中 是维已知向量 是维解向量 当时 即方程组中方程的个数多于未知量的个数时 称此方程组为超定方程组 2 2线性方程组求解 一般来说 超定方程组无解 此时为矛盾方程组 这时需要寻找方程组的一个 最近似 的解 如最小二乘解 线性方程组 2 2 可能无解 即任何一组数 都可能使不等于零 我们设法找 使 2 3 最小 称为方程组的最小二乘解 2 2线性方程组求解 最小二乘法 最小二乘解所满足的代数方程 又称为正规方程组 为 2 4 且线性方程组 2 4 总是有解 2 2线性方程组求解 例2 9求超定方程组的最小二乘解 解 原方程组写成矩阵形式为 则正规方程组为 2 2线性方程组求解 即令 利用MATLAB中矩阵的左除 X A b即得 2 2线性方程组求解 例2 10已知图2 1给出了某城市部分单行街道的交通流量 每小时通过的车辆数 图中共有6个路口 已有9条街道记录了当天的平均车流量 另有7条街道的平均车流量未知 试根据每个路口的进出车流量相等关系推算这7处的平均车流量 图2 1 2 3应用实例 解 在每一个路口处可根据进出的车流量相等关系 建立一个线性代数方程 由图2 1中有六个路口 可建立含六个方程的线性方程组 将方程组写成矩阵向量形式为 2 3应用实例 求解问题分为三个步骤 第一步 判断方程组是否有解 第二步 如果有解则求出方程组的通解 第三步 在通解中找非负特解 实际问题所需 A 1010000 1 101000 0100 100 0010010 000101 1 000010 1 b 700 200 200 500 0 200 2 3应用实例 r1 rank A r2 rank Ab r1 5r2 5 rref Ab 2 3应用实例 求系数矩阵和增广矩阵的秩 Ab 10000 10200010000 100010010500000101 10000010 1 20000000000 2 3应用实例 2 3应用实例 由最简行阶梯形矩阵 可得与原方程组等价的简化后的方程组 取为自由未知数 可得原方程组的通解为 其中 为任意实数 对方程组通解表达式 取适当的和 使特解为非负数 即得一组满足问题条件的解 2 3应用实例 如 取 得 2 3应用实例 即求出路口7处的平均车流量的一个解 例2 12 世界人口预测问题 据统计 六十年代世界人口数据如下 有人根据表中数据 预测公元2008年世界人口会超过70亿 作出这一预测结果所用的方法就是数据拟合方法 根据数学模型 构造出能逼近表2 2中数据的拟合函数 根据表中的数据来构造拟合函数 预测公元2008年时的世界人口数 2 3应用实例 解 据人口增长的统计资料和人口理论数学模型 当人口总数N不是很大时 在不太长的时期内 人口增长接近于指数增长 因此 采用指数函数 接下来对数据进行拟合 将上式两边同取对数 得 2 3应用实例 拟合问题要先确定函数模型 插值问题则不确定 令或 拟合函数则变换 由计算 得世界人口的拟合函数的人口数与年份的对应关系表 2 3应用实例 根据表中数据及等式 可列出关于两个未知数的9个方程的超定方程组 2 3应用实例 输入 t 1960 1968输出 t 196019611962196319641965196619671968 2 3应用实例 MATLAB求解过程如下 y 3 39183 42133 45033 46983 47633 49203 51333 53223 5505 A ones 9 1 t 2 3应用实例 输出超定方程的系数矩阵 A 119601196111962119631196411965119661196711968 2 3应用实例 A A A y ab 33 04310 0186即得a 33 0383 b 0 0186错的 将它代入拟合函数 2 3应用实例 由A A ab A y 输入 exp 33 0383 0 0186 2008 输出 N 2008 74 4777 即2008年的世界人口预测为74 4777亿 2 3应用实例 图形输出 输入 t 1960 1968 t0 2008 N 29 7230 6131 5132 1332 3432 8533 5634 2034 83 y log N A ones 9 1 t d A y a d 1 b d 2 N0 exp a b t0 x 1960 2010 yy exp a b x plot x yy t N o 2008 N0 o 2 3应用实例 2 3应用实例 输出人口曲线 1 方程求根2 线性方程组求解3 插值 拟合4 应用 人口问题 划艇比赛成绩 主要内容 3 插值与拟合 2 拟合的基本原理 线性最小二乘拟合 3 面对实际问题 应该用插值还是拟合 2020年3月21日 67 插值问题实例1 标准正态分布函数 x 求 1 014 1 014 0 8438 0 8461 0 8438 0 4 0 8447 68 插值问题实例2 机械加工 69 插值问题的提法 70 求解插值问题的基本思路 71 1 拉格朗日 Lagrange 仅用于多项式插值 1 1插值多项式 有唯一解 72 1 拉格朗日 Lagrange 多项式插值 1 2拉格朗插值多项式 又 2 有唯一解 故 3 与 1 相同 73 1 拉格朗日 Lagrange 多项式插值 1 3误差估计 74 1 拉格朗日 Lagrange 多项式插值 1 4例将 0 2 n等分 用g x cos x 产生n 1个节点 作Ln x 取n 1 2 计算cos 6 估计误差 解 n 1 x0 y0 0 1 x1 y1 2 0 L1 x y0l0 y1l1 1 2x cos 6 0 6667 n 2 x0 y0 0 1 x1 y1 4 0 7071 x2 y2 2 0 L2 x y0l0 y1l1 y2l2 8 x 4 x 2 2 16x x 2 0 7071 2cos 6 L2 6 0 8508精确值 cos 6 0 8660 75 1 拉格朗日 Lagrange 多项式插值 1 5一般不用拉格朗日插值多项式的振荡 Runge 融合现象 76 2 分段线性插值 得出的曲线不连续 很致命 计算量与n无关 n越大 误差越小 77 3 三次样条插值 插值的MATLAB命令 interp1 1 D插值 interp2 2 D插值 interp3 3 D插值 interpn N D插值 Syntax yi interp1 x Y xi yi interp1 Y xi yi interp1 x Y xi method yi interp1 x Y xi method extrap yi interp1 x Y xi method extrapval pp interp1 x Y method pp 插值的MATLAB命令 X Y meshgrid 3 25 3 Z peaks X Y XI YI meshgrid 3 125 3 ZI interp2 X Y Z XI YI mesh X Y Z hold mesh XI YI ZI 15 holdoffaxis 33 33 520 无曲线类型拟合时的插值实例 某观测站测得某日6 00至18 00之间每隔2小时的室内外温度 用线性插值求得该日室内外6 30至17 30之间每隔2小时各点的近似温度 h 6 2 18 t 18 20 22 25 30 28 24 15 19 24 28 34 32 30 XI 6 5 2 17 5 YI interp1 hX tY XI linear 线性 YI 18 500016 000020 500020 250022 750025 000026 250029 500029 500033 500027 000031 5000 插值举例 经过插值达到数据更为平滑的效果举例 有一组海底高程测量数据 采用分段线性插值方式绘制海底形状图 因为实测数据有限 所以绘制的海底地形图很不光滑 如图所示 surf X Y Z view 25 25 绘图 视角 X Y为平面坐标 Z为实测高程数据 插值举例 为了进行插值 将计算区域由网格生成函数meshgrid生成更密的网格 经过插值后海底地形图及代码如下 x linspace 5 5 50 y linspace 5 5 50 XI YI meshgrid x y 在X Y坐标范围内生成更密网格 ZI interp2 X Y Z XI YI cubic 曲线插值 surf XI YI ZI view 25 25 绘图 视角 插值举例 2020年3月21日 83 曲线拟合问题的提法 已知一组 二维 数据 即平面上n个点 xi yi i 1 n 寻求一个函数 曲线 y f x 使f x 在某种准则下与所有数据点最为接近 即曲线拟合得最好 y f x i为点 xi yi 与曲线y f x 的距离 2020年3月21日 84 曲线拟合问题最常用的解法 线性最小二乘法的基本思路 先选定一组函数r1 x r2 x rm x m n 令f x a1r1 x a2r2 x amrm x 1 其中a1 a2 am为待定系数 确定a1 a2 am的准则 最小二乘准则 使n个点 xi yi 与曲线y f x 的距离 i的平方和最小 记 问题归结为 求a1 a2 am使J a1 a2 am 最小 2020年3月21日 85 线性最小二乘法的求解 记 当RTR可逆时 4 有唯一解 2020年3月21日 86 线性最小二乘拟合f x a1r1 x amrm x 中函数 r1 x rm x 的选取 1 通过机理分析建立数学模型来确定f x 2 将数据 xi yi i 1 n作图 通过直观判断确定f x 87 用MATLAB作线性最小二乘拟合 仅能对多项式拟合 1 作多项式f x a1xm amx am 1拟合 可利用已有程序 a polyfit x y m a是一个向量 输入 数据x y 同长度数组 m 拟合多项式次数 输出 系数a a1 am am 1 数组 2 对超定方程组 可得最小二乘意义下的解 仍用 拟合 多项式曲线拟合用polyfit函数拟合数据集的多项式在最小二乘意义上的系数 调用格式为 P polyfit x y n x和y是包含要拟合的x和y数据的矢量 n是多项式的阶次 例 x 12345 y 5 543 1128290 7498 4 p polyfit x y 3 p 0 191731 5821 60 326235 3400 拟合 x2 1 1 5 y2 polyval p x2 polyval 求值 plot x y o x2 y2 gridon 拟合 91 用MATLAB作线性最小二乘拟合 已知温度和热敏电阻的实验数据如下表所示 并已知二者成直线关系 求t 60度时电阻R t 20 5 32 7 51 73 95 7 R 765 826 873 942 1032 p polyfit t R 1 1次多项式 y ax b 其中a p 1 b p 2 p 3 3987393392826702 096806861956 x min t 0 1 max t y polyval p x 计算多项式的值 plot x y b t R or 绘图 R60 polyval p 60 计算60度时R的值R60R60 906 021167218912在这种情况下 计算得到的值比较准确 可以消除观测带来的误差 而如果改为插值 如果用线性插值方法 则只与相邻两个实测数据有关 对计算结果的误差影响较大 用MATLAB作线性最小二乘拟合 1 方程求根2 线性方程组求解3 插值 拟合 最小二乘仅解决线性问题 4 应用 人口问题 划艇比赛成绩 主要内容 建模应用 1 人口问题已知美国1790 1900年的人口数据 估计指数模型 中的参数N0与r N0与r不是线性关系 对指数模型两边取对数 得 记y lnN t a lnN0 则有 a与r是线性关系 程序如下 x 0 11 离散的时间点y 3 9 5 9 7 2 9 6 12 9 17 1 23 2 31 4 38 6 50 2 62 9 76 对应时间点的数据z log y a polyfit x z 1 线性多项式插值思考 如果你拟合的数据一直取到2000年 得到还是这个拟合方程么 2 划艇比赛通过课上的比例分析我们推得 题中所给的表格赛艇2000米成绩t 分 种类1234平均单人7 167 257 287 177 21双人6 876 926 956 776 88四人6 336 426 486 136 32八人5 875 925 825 735 84 利用4次国际大赛冠军的平均成绩对模型tn 1 9进行检验 与模型巧合 程序 y 7 216 886 325 84 x 1248 z log x t log y a poly z t 1 用MATLAB作插值计算 已知下轮廓线上数据如下 求x每改变0 1时的y值 上机练习 2020年3月21日 101 求600C时的电阻R 设R at ba b为待定系数 上机练习 2020年3月21日 102 求血药浓度随时间的变化规律c t 半对数坐标系 semilogy 下的图形 上机练习

展开阅读全文