用Matlab进行现代科学运算.ppt

Matlab基础与应用 第三讲 2 用Matlab进行现代科学计算 张文静zhangwj 解析解与数值解数值线性代数问题及求解数值微积分问题常微分方程的数值解法 主要内容 解析解与数值解数值线性代数问题及求解数值微积分问题常微分方程的数值解法 主要内容 5 1解析解与数值解 在实际的工程技术中 一般通过数值解法来获得问题的解 至少有两个原因 1 解析解不存在如圆周率的值就没有解析解 超越积分也没有解析解 在这种情况下 就必须采用数值解技术 2 解析解存在但不实用如求解n元一次代数方程组的问题 从理论上讲 总可以把多元一次方程组简化成解析可解的形式 然而当n较大时 需要的基本运算次数是非常惊人的 对于求解实际问题所需要的计算成本是根本无法接受的 只能采用数值解法 5 1解析解与数值解 解析解与数值解数值线性代数问题及求解数值微积分问题常微分方程的数值解法 主要内容 5 2 1特殊矩阵的MATLAB输入 零矩阵 zeros m n m n 定义矩阵维数 幺矩阵 ones m n 单位矩阵 eye m n 随机元素矩阵 rand m n 0 1 上均匀分布randn m n 正态分布 对角阵 diag v v为对角元素组成的向量 5 2 2矩阵的特征参数运算 矩阵的行列式 det A 矩阵的迹 trace A 即对角线元素之和 矩阵的秩 rank A 即线性无关的列数或行数 矩阵的范数norm A p 当A为向量时 范数定义稍有不同 5 2 2矩阵的特征参数运算 矩阵的特征多项式 c poly A 矩阵的特征值 特征方程的根roots poly A 或eig A 若A为矩阵 则返回的c就是向量 1 c1 cn 1 cn 降幂排列 5 2 2矩阵的特征参数运算 矩阵的特征值与特征向量 v d eig A 其中 矩阵v的各列为特征向量 矩阵d的对角元素为特征值 a 123 456 780 v d eig a v 0 2998 0 7471 0 2763 0 70750 6582 0 3884 0 6400 0 09310 8791d 12 1229000 0 3884000 5 7345 5 2 2矩阵的特征参数运算 矩阵求逆 B inv A 为求矩阵A的逆矩阵 a 123 456 780 b inv a c a b d b a bcd ans 1 77780 8889 0 11111 00000 0 00001 00000 000001 5556 0 77780 2222 0 00001 00000 0 00001 00000 0 11110 2222 0 11110 0000 0 00001 00000 00000 00001 000 5 2 2矩阵的特征参数运算 多项式及多项式矩阵求值 多项式求值B polyval aa x 多项式矩阵求值B polyvalm aa A 其中aa为多项式系数降幂排列构成的向量 x为指定的标量 A为给定的方阵 5 2 2矩阵的特征参数运算 面向矩阵各个元素的函数 B 函数名 A 用命令helpelfun可以查看这些命令列表 主要有 sinsinhasinasinhcoscoshacosacoshtantanhatanatan2atanhsecsechasecasechcsccschacscacschcotcothacotacothexploglog10log2pow2sqrtabsanglecomplexconjimagrealunwrapisrealcplxpairfixfloorceilroundmodremsign 5 2 2矩阵的特征参数运算 对矩阵进行数值分析的函数 B 函数名 A 用命令helpdatafun可以查看这些命令列表 主要有 maxminmeanmedianstdvarsortsortrowssumprodhisthistctrapzcumsumcumprodcumtrapzdiffgradientdel2corrcoefcovsubspacefilterfilter2convconv2convndeconvdetrendfftfft2fftnifftifft2ifftnfftshiftifftshift 5 2 2矩阵的特征参数运算 5 2 3矩阵的相似变换与分解 1 三角形分解最基本的分解方法 可把任何一个方阵表示成两个三角矩阵的乘积 其中一个方阵是换位的下三角矩阵 另一个是上三角矩阵 这种分解法常称为LU分解或LR分解 其算法多为高斯消去法 是求解线性方程组 inv 和det 的基础 用lu 函数可得到分解后的两个三角矩阵 其调用形式为 L U lu A 2 正交分解也称为 QR 分解 将一个矩阵表示成一个正交矩阵和一个上三角矩阵的乘积 A Q R用qr 函数可得到分解后的两个矩阵 其调用形式为 Q R qr A 5 2 3矩阵的相似变换与分解 3 奇异值分解 U S V svd A 其中 U和V是正交矩阵 而S是对角矩阵 即A U S V 奇异值分解是矩阵分析的有力工具 Matlab的一些函数也是基于此实现的 包括求伪逆函数pinv A 求秩函数rank A 求欧式范数norm A 2 求条件数cond A 5 2 3矩阵的相似变换与分解 解析解与数值解数值线性代数问题及求解数值微积分问题常微分方程的数值解法 主要内容 MATLAB提供了计算给定向量差分的函数diff 其调用方法是dy diff y 按列运算 行数减1 5 3 1数值差分运算 diff y magic 6 y 351626192433272123253192222720828331710153053412141643629131811 diff y ans 32311 54128 23 514 5 231931 5 17 522 231 541 2631 514 5 5 3 1数值差分运算 5 3 2数值积分 求解函数定积分的数值方法是多种多样的 如梯度法 Simpson法 Romberg法等都是经常采用的方法 基本思想 将整个积分空间分割成若干个自空间 而每个小空间上的函数积分可求 因而整个空间函数积分可求 y n sum F a b 等宽矩形法求定积分 y n trapz F a b 梯形法求定积分 y n quad F a b tol 自适应辛普森积分 y n quadl F a b tol 自适应Lobatto积分 牛顿8段积分 替代quad8 其中F为描述被积函数的字符串变量 一般为一个F m函数文件名 加引号 还可用inline 来定义一个函数 a b为积分上下限 tol为变步长积分用误差限 y为积分结果 n为被积分函数的调用次数 常见的一元函数数值积分指令 5 3 2数值积分 例 求无穷定积分 f inline 1 sqrt 2 pi exp x 2 2 x y kk quad f 8 8 y 1 00000197533430kk 81 y1 kk1 quadl f 8 8 y1 1 00000000000023kk1 161 y kk quad f 15 15 y 0 99999920879563kk 89 y1 kk1 quadl f 15 15 y1 0 99999999999999kk1 769 该无穷定积分的理论值为1 5 3 2数值积分 例子 f inline exp x 2 2 sin x 2 y 2 pi x y y dblquad f 2 2 1 1 y 0 2506 双重积分函数 dblquad F x m x M y m y M tol 三重积分函数 triplequad F x m x M y m y M z m z M tol 5 3 2数值积分 解析解与数值解数值线性代数问题及求解数值微积分问题常微分方程的数值解法 主要内容 通常把含有自变量t 未知的一元函数x t 及其导数或微分的方程 叫做常微分方程 ODE OrdinaryDifferentialEquation 如一阶常微分方程的初值问题 5 4 1一般常微分方程的数值解法 求解常微分方程组的数值方法是多种多样的 如常用的Euler法 Runge Kutta法 Adams线性多步法 Gear法等 MATLAB求解常微分方程的函数如下 5 4 1一般常微分方程的数值解法 ode23 和ode45 最常用 采用自适应变步长求解方法 t x ode23 方程函数名 t0 tf x0 选项 t x ode45 方程函数名 t0 tf x0 选项 方程函数名为描述系统状态方程的M函数的名称 该函数名应该用引号括起来 t0和tf分别为起始和终止计算时间 x0为系统的初始状态变量的值 返回值为求解的时间变量和相应的状态列向量构成的矩阵转置 5 4 1一般常微分方程的数值解法 方程函数名的编写格式是固定的 functionxdot 方程函数名 t x 其中t为时间变量 x为方程的状态变量 xdot代表状态变量的导数 注意 即时微分方程是非时变的 也应该在函数输入变量列表中写上t占位 5 4 1一般常微分方程的数值解法 5 4 1一般常微分方程的数值解法 例 解Lorenz模型状态方程 在lorenzeq m函数中 functionxdot lorenzeq t x xdot 8 3 x 1 x 2 x 3 10 x 2 10 x 3 x 1 x 2 28 x 2 x 3 t final 100 x0 0 0 1e 10 t x ode45 lorenzeq 0 t final x0 figu