MATLAB第5章数值计算.ppt

《MATLAB第5章数值计算.ppt》由会员分享，可在线阅读，更多相关《MATLAB第5章数值计算.ppt（43页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、第五章第五章 数值计算数值计算主要内容主要内容5.1 线性代数线性代数5.2 函数分析函数分析5.3 数据拟合数据拟合5.4 插值和样条插值和样条5.5 常微分方程的数值解常微分方程的数值解5.1 线性代数线性代数 一、一、LU分解分解 1.LU分解分解 一个矩阵可以分解为一个上三角矩阵和一个下三角矩阵的乘积，一个矩阵可以分解为一个上三角矩阵和一个下三角矩阵的乘积，称之为称之为LU分解。分解。LU分解是用高斯主元消去法实现的，通常要对主元分解是用高斯主元消去法实现的，通常要对主元位置进行交换，主元交换的方法是将被分解矩阵左乘一个由位置进行交换，主元交换的方法是将被分解矩阵左乘一个由01构成构成

2、的行交换阵。的行交换阵。【调用格式调用格式】L,U,P=lu(X)对矩阵对矩阵X进行进行LU分解，并进行主元交换分解，并进行主元交换L,U=lu(X)对矩阵对矩阵X进行进行LU分解，无主元交换分解，无主元交换 【说明说明】X为被分解的矩阵，为被分解的矩阵，L为主对角元素为为主对角元素为1的下三角矩的下三角矩 阵，阵，U为上三角矩阵，为上三角矩阵，P为行交换矩阵。为行交换矩阵。5.1 线性代数线性代数2.行列式和求逆行列式和求逆矩矩阵的行列式和求逆可以通的行列式和求逆可以通过LU分解的方法求解，分解的方法求解，Matlab提供了相关函数。提供了相关函数。【调用格式用格式】d=det(X)求矩求矩

3、阵X的行列式的行列式 Y=inv(X)求矩求矩阵X的逆矩的逆矩阵例例例例5.1.15.1.15.1 线性代数线性代数 二、特征二、特征值和特征向量和特征向量 对于求解矩于求解矩阵的特征的特征值和特征向量，和特征向量，Matlab提供了提供了eig函数。函数。【调用格式用格式】D=eig(A)计算矩算矩阵A的特征的特征值，d为特征特征值构成的向量构成的向量 V,D=eig(A)计算矩算矩阵A的特征的特征值对角角阵D和特征向量矩和特征向量矩阵V V,D=eig(A,nobalance)当矩当矩阵A中有与截断中有与截断误差近似的数差近似的数值，用本指令，用本指令例例例例5.1.25.1.2三、奇异三

4、、奇异值分解分解1.矩阵的奇异值分解矩阵的奇异值分解任意矩阵任意矩阵A可进行奇异值分解，即存在酉矩阵可进行奇异值分解，即存在酉矩阵U和和V，使下面等式成立，使下面等式成立其中其中 称为矩阵的奇异值。称为矩阵的奇异值。5.1 线性代数线性代数【调用格式调用格式】s=svd(A)求矩阵求矩阵A的奇异值，的奇异值，s为由奇异值为由奇异值 构成的向量构成的向量U,S,V=svd(A)矩阵矩阵A的奇异值分解的奇异值分解5.1 线性代数线性代数2.矩阵结构特征的奇异值描述矩阵结构特征的奇异值描述 矩阵的奇异值可以描述矩阵的结构特征。有关矩阵结构特征的矩阵的奇异值可以描述矩阵的结构特征。有关矩阵结构特征的M

5、ATLAB 函数有如下几种。函数有如下几种。r=rank(A,tol)在指定容差在指定容差tol下，求矩阵下，求矩阵A的秩。的秩。tol可以省略。可以省略。Z=null(A)求矩阵求矩阵A的零空间。的零空间。V=orth(A)求矩阵求矩阵A的值空间。的值空间。n=norm(A)求矩阵求矩阵A的的2范数。范数。n=norm(A,p)求矩阵求矩阵A的各种范数。的各种范数。c=cond(X,p)求矩阵求矩阵A的条件数，的条件数，p可以省略。可以省略。theta=subspace(A,B)求求A和和B矩阵所张子空间的夹角。矩阵所张子空间的夹角。B=pinv(A,tol)在指定容差在指定容差tol下，求

6、矩阵下，求矩阵A的的 广义逆，广义逆，tol可以省略。可以省略。5.1 线性代数线性代数 四、四、线性方程性方程组的解的解形如形如 的的线性方程性方程组中中独立方程的个数等于独立未知参数的个数称独立方程的个数等于独立未知参数的个数称为恰定方程；恰定方程；独立方程的个数大于独立未知参数的个数称独立方程的个数大于独立未知参数的个数称为超定方程；超定方程；独立方程的个数小于独立未知参数的个数称独立方程的个数小于独立未知参数的个数称为欠定方程。欠定方程。针对不同情况，可以采用以下不同情况，可以采用以下3种方法求解：种方法求解：1.左除运算符法左除运算符法 对于一般的非奇异矩于一般的非奇异矩阵A，可以求

7、得唯一数，可以求得唯一数值解。欠定解。欠定 方程和超定方程，可以方程和超定方程，可以获得最小二乘解。得最小二乘解。x=Ab5.1 线性代数线性代数2.广广义逆法逆法如果用左除运算符求解的如果用左除运算符求解的时候出候出现提示矩提示矩阵A为非奇异的警告或者解中非奇异的警告或者解中出出现Nan，则可以采用广可以采用广义逆法。逆法。x=pinv(A)*b3.符号符号计算法算法 可以求得方程可以求得方程组的符号解，的符号解，对于欠定方程可以求得具有自由于欠定方程可以求得具有自由变量的解。量的解。5.1 线性代数线性代数例例5.1.3 求以下求以下3个方程个方程组的解的解 I：II：III：5.2 函数

8、分析函数分析 一、函数的零点一、函数的零点1.多多项式的根式的根 通通过roots函数来求取多函数来求取多项式全部的根。式全部的根。【调用格式用格式】r=roots(p)多多项式求根函数，求取多式求根函数，求取多项式的全部式的全部根根【说明明】p为多多项式的系数行向量，式的系数行向量，r为多多项式所有根构成的列向量式所有根构成的列向量5.2 函数分析函数分析2.一元函数零点一元函数零点【调用格式调用格式】x,fval,exitflag,output=fzero(fun,x0,options)完整调用格式完整调用格式 x=fzero(fun,x0)一元函数零点的最简调用格式一元函数零点的最简调用

9、格式 【说明说明】fzero只能求得只能求得x0附近的单个零点，不能求取函数附近的单个零点，不能求取函数的所有零点。的所有零点。输入变量输入变量fun表示一元函数，可以是字符串、内联表示一元函数，可以是字符串、内联函数或者函数句柄。函数或者函数句柄。输入变量输入变量x0为零点的初始猜测值（自变量值）。为零点的初始猜测值（自变量值）。如果如果x0为标量，则求距离为标量，则求距离x0最近的那个零点；如最近的那个零点；如果果x0=a,b，要求，要求fun(a)和和fun(b)异号，此时求自异号，此时求自变量在变量在a,b区间内的零点。区间内的零点。5.2 函数分析函数分析输入变量输入变量option

10、s为优化迭代选项，是一个结构体。为优化迭代选项，是一个结构体。输出变量输出变量x为零点处的自变量值，输出变量为零点处的自变量值，输出变量fval为为 零点处的函数值。零点处的函数值。输出变量输出变量exitflag表示函数中止计算的条件。表示函数中止计算的条件。若若exitflag0表示找到零点后退出。表示找到零点后退出。例例5.2.1 求函数求函数在在区区间间内的所有零点。内的所有零点。3.多元函数的零点多元函数的零点多元函数的精确零点可以用多元函数的精确零点可以用fsolve函数求取，函数求取，但是必须提供零点的大致位置才能进行数值搜索。但是必须提供零点的大致位置才能进行数值搜索。5.2

11、函数分析函数分析【调用格式调用格式】x,fval,exitflag,output=fsolve(fun,x0,options)完整格式完整格式x=fsolve(fun,x0)多元非线性方程求解的最简格式多元非线性方程求解的最简格式【说明说明】fun表示自变量为向量的函数，可以是字符串、表示自变量为向量的函数，可以是字符串、内联函数、函数句柄。内联函数、函数句柄。输入变量输入变量x0为零点的初始猜测向量（自变量值）。为零点的初始猜测向量（自变量值）。输出变量输出变量x为零点处的自变量值，输出变量为零点处的自变量值，输出变量fval为零为零 点处的函数值。点处的函数值。例例5.2.2 求方程求方程

12、组组在在附近的解。附近的解。5.2 函数分析函数分析 二、函数的极值点二、函数的极值点Matlab提供了提供了3个求极值点的函数，其输入输出参数的定义和个求极值点的函数，其输入输出参数的定义和 fsolve函数基本相同。函数基本相同。1.求一元函数求一元函数fun在自变量（在自变量（x1，x2）区间的最小值）区间的最小值 x,fval,exitflag,output=fminbnd(fun,x1,x2,options)2.用单纯形法求多元函数用单纯形法求多元函数fun在自变量向量在自变量向量x0附近的极小值点附近的极小值点 x,fval,exitflag,output=fminsearch(f

13、un,x0,options)3.用拟牛顿法求多元函数用拟牛顿法求多元函数fun在自变量向量在自变量向量x0附近的极小值点附近的极小值点 x,fval,exitflag,output=fminunc(fun,x0,options)例例5.2.3 求二元函数求二元函数在在附近的极小附近的极小值值。5.2 函数分析函数分析三、数三、数值微分微分数数值导数和微分是基于差分定数和微分是基于差分定义的，的，对原始数据的采原始数据的采样间隔依隔依赖很大，很大，如果原始数据含有噪声，如果原始数据含有噪声，则数数值导数数结果没有什么价果没有什么价值，因此尽量避免求，因此尽量避免求数数值导数。数。Matlab中用

14、中用diff函数求相函数求相邻数据数据间的数的数值差分。差分。【调用格式用格式】DX=diff(X)求求X相相邻元素的一元素的一阶差分，即后一个元素减去当前元素差分，即后一个元素减去当前元素DX=diff(X,n)求求X相相邻元素的元素的n阶差分差分DX=diff(X,n,dim)在在dim指定的指定的维上，求上，求X 相相邻元素的元素的n阶差分差分5.2 函数分析函数分析【说明说明】如果如果X是向量，差分是相邻元素相减；如果是向量，差分是相邻元素相减；如果X是矩阵，差分是相邻是矩阵，差分是相邻 行行 间对应元素相减。间对应元素相减。差分数据差分数据DX元素的个数要比被操作数元素的个数要比被操

15、作数X少一个。少一个。DX只是差分数据，如果相邻数据点之间的步长为只是差分数据，如果相邻数据点之间的步长为DH，则，则DX./DH为为 导数数据。导数数据。例例5.2.4 绘绘制函数制函数的的导导函数在函数在区区间间的曲的曲线线。5.2 函数分析函数分析四、数四、数值积分分数数值积分分分分为开型开型积分和分和闭型型积分，二者的区分，二者的区别在于是否在于是否计算算积分分区区间端点端点处的函数的函数值。Matlab提供的数提供的数值积分函数有些适用于开型分函数有些适用于开型积分分和和闭型型积分，有些只能用于分，有些只能用于闭型型积分。分。积分分计算可以采用本算可以采用本节介介绍的相的相关函数，也

16、可以采用关函数，也可以采用样条条积分法，分法，还可以采用符号可以采用符号积分方法。分方法。1.一元函数的数一元函数的数值积分分求一元函数数求一元函数数值积分的函数有很多，每个函数采用不同的数分的函数有很多，每个函数采用不同的数值算法，算法，各有各有优缺点，精度和运算速度也不尽相同。缺点，精度和运算速度也不尽相同。本本节只介只介绍最常用的数最常用的数值积分函数分函数quadl。【调用格式调用格式】q=quadl(fun,a,b,tol,trace)采用递推自适应采用递推自适应Lobatto法计算一元函数的积分法计算一元函数的积分5.2 函数分析函数分析【说明说明】fun为被积函数，可以用字符串、

17、内联函数和函数为被积函数，可以用字符串、内联函数和函数M文件的函数文件的函数 句柄表示，且被积函数表达式要按照数组运算规则来编写。通常句柄表示，且被积函数表达式要按照数组运算规则来编写。通常积分变量采用字母积分变量采用字母x。a和和b为积分变量的积分上下限，为常数数值。为积分变量的积分上下限，为常数数值。tol为绝对误差，是一个标量，可以省略。为绝对误差，是一个标量，可以省略。trace为跟踪标志，当为跟踪标志，当trace为非零值时，随积分进程会为非零值时，随积分进程会 逐点画出被积函数。逐点画出被积函数。返回值返回值q为数值积分的结果。为数值积分的结果。例例例例5.2.5 5.2.5 求求

18、求求积积积积分分分分形如形如 的的闭闭型二重数型二重数值积值积分的函数分的函数为为dblquad【调调用格式用格式】q=dblquad(fun,x1,x2,y1,y2,tol,method)2.二重数二重数值积值积分分【说说明明】fun为为被被积积函数，函数，可以用字符串、内可以用字符串、内联联函数和函数函数和函数M文件的函数句柄文件的函数句柄 表示。表示。积积分分变变量用量用x，y表示，表示，x为为内重内重积积分分变变量，量，y为为外重外重积积分分变变量。量。x2和和x1是是变变量量x的的积积分上下限，分上下限，y2和和y1是是变变量量y的的积积分上下限。分上下限。method表示表示选选用

19、的用的积积分算法，可以省略，确是算法的分算法，可以省略，确是算法的 是是quad，还还可以可以选则选则quadl或者用或者用户户自定自定义义的的积积分分 算法算法M函数文件的函数句柄。函数文件的函数句柄。本函数无法本函数无法计计算内重算内重积积分上下限分上下限为为函数表达式的情况。函数表达式的情况。5.2 函数分析函数分析求形如求形如 的的闭闭型三重数型三重数值积值积分的函数分的函数为为triplequad。【调调用格式用格式】q=triplequad(fun,x1,x2,y1,y2,z1,z2,tol,method)3.三重数三重数值积值积分分【说说明明】fun为为被被积积函数函数 ，可以用

20、字符串、内，可以用字符串、内联联函数和函数函数和函数M文件的函数句文件的函数句 柄表示。柄表示。积积分分变变量用量用x，y，z表示，表示，从内到外从内到外积积分分变变量依次量依次为为x，y，z。x2和和x1是是变变量量x的的积积分上下限，分上下限，y2和和y1是是变变量量y的的积积分分 上下限，上下限，z2和和z1是是变变量量z的的积积分上下限。分上下限。其他事其他事项项同函数同函数dblquad。5.2 函数分析函数分析例例5.2.6 计计算定算定积积分分5.3 数据拟合数据拟合已知某些离散的原始数据，建立一个曲已知某些离散的原始数据，建立一个曲线方程，方程，让它以最佳的方它以最佳的方式反映

21、原始数据的式反映原始数据的变化化趋势，尽量避免出，尽量避免出现局部的波局部的波动，这种方法称种方法称为拟合。不要求合。不要求拟合曲合曲线经过所有的原始数据。最佳方式通常是指用所有的原始数据。最佳方式通常是指用拟合得到的数学模型合得到的数学模型计算出来的算出来的计算数据和原始数据之算数据和原始数据之间的的误差的平差的平方和最小，方和最小，这种方式称种方式称为最小二乘。通常是已知曲最小二乘。通常是已知曲线方程的形式（数方程的形式（数学模型的学模型的结构），根据原始数据来构），根据原始数据来计算曲算曲线方程的参数（数学模型的方程的参数（数学模型的参数）。本参数）。本节介介绍Matlab中如何中如何实

22、现最小二乘最小二乘拟合。合。5.3 数据拟合数据拟合一、多一、多项式式拟合合多多项式式拟合是用一个合是用一个n阶多多项式模型去式模型去拟合原始数据的方法，需要合原始数据的方法，需要计算的算的模型参数包括多模型参数包括多项式的式的阶次和多次和多项式的系数。式的系数。Matlab提供了提供了polyfit函数函数实现多多项式式拟合。合。【调用格式用格式】p=polyfit(x,y,n)根据根据给定的数据定的数据(x,y)，计算算n阶拟合多合多项式的系数向量式的系数向量p ye=polyval(p,x)计算自算自变量量为x时多多项式式p的的值（估（估计值）5.3 数据拟合数据拟合 【说明说明】多项式

23、拟合一般不要超过多项式拟合一般不要超过5阶，否则计算误差变大。阶，否则计算误差变大。拟合多项式只在原始数据范围内可以保证精度，超出范围使用拟合多项式只在原始数据范围内可以保证精度，超出范围使用 拟合多项式无法保证预报的精度。拟合多项式无法保证预报的精度。例例5.3.1：已知五个数据点：已知五个数据点：1,5.5,2,43,3,128,4,290,5,498，试画出这五个点拟合的三次曲线。试画出这五个点拟合的三次曲线。二、最小二乘二、最小二乘拟合合1.线性最小二乘性最小二乘拟合合对于于线性数学模型的参数估性数学模型的参数估计，可以用形如，可以用形如Y=Ax+b的一的一阶多多项式式拟合来估合来估计

24、参数。某些非参数。某些非线性模型性模型经过变量替量替换也可以也可以转换为线性模型，也性模型，也可以采用可以采用线性估性估计方法。方法。【调用格式用格式】a=lsqlin(Fun,a0,lb,ub,options,p1,p2,.)a=lsqnonlin(Fun,a0)5.3 数据拟合数据拟合例例5.3.2：对于非线性数学模型：对于非线性数学模型 ，其中，其中a和和b为模型为模型参数，参数，采用变量替换将其转换为线性数学模型。采用变量替换将其转换为线性数学模型。解：对已知数学模型取自然对数有解：对已知数学模型取自然对数有 令令 则有的线性数学模型则有的线性数学模型5.3 数据拟合数据拟合函数的功能

25、是求取向量的函数的功能是求取向量的值值，使向量函数，使向量函数 满足满足输入参数输入参数Fun是被解函数，是被解函数，Fun是向量函数，其自变是向量函数，其自变量量a是向量。是向量。Fun可以用字符串、内联函数或者函数可以用字符串、内联函数或者函数 M文件的函数句柄表示。文件的函数句柄表示。2.非非线线性最小二乘估性最小二乘估计计可以使用可以使用MATLAB提供的提供的lsqnonlin函数函数实现实现非非线线性最小二乘估性最小二乘估计计。【调调用格式用格式】a,resnorm,residual,exitflag,output=lsqnonlin(Fun,a0,lb,ub,options,p1

26、,p2,.)a=lsqnonlin(Fun,a0)【说说明明】5.3 数据拟合数据拟合输入变量输入变量a0为向量，表示所求变量为向量，表示所求变量a的初始猜测值。的初始猜测值。输输入入变变量量lb表示所求表示所求变变量量a的下限，的下限，输输入入变变量量ub表示所求表示所求变变量量a的上的上限。如果没有上下限限制，可以用空矩限。如果没有上下限限制，可以用空矩阵阵表示，如果没有下限表示，如果没有下限可以用可以用-inf表示。表示。从第六个从第六个输输入入变变量开始的参数量开始的参数p1,p2等是向等是向Fun函数函数传递传递的参数，的参数，其名字其名字应该应该和定和定义义Fun函数的函数的输输入

27、入变变量名相一致。量名相一致。输输出出变变量量a表示所求表示所求变变量的量的值值。如果将如果将 构造构造为为残差函数，且残差函数，且 为为被估被估计计参数参数向量，向量，则则lsqnonlin就就变为变为求解非求解非线线性最小二乘估性最小二乘估计计的函数。的函数。5.3 数据拟合数据拟合例例5.3.3 混凝土的抗混凝土的抗压压强强度度 随养随养护时间护时间 的延的延长长而增加，其数学模型而增加，其数学模型为为现现将一批混凝土作成将一批混凝土作成12个个测试块测试块，记录记录了养了养护时间护时间x（日）及抗（日）及抗压压强强度度y （kg/cm2）的数据，如表）的数据，如表5.3.1 所示。所示

28、。试试估估计该计该数学模型的参数数学模型的参数a1，a2，a3，a4。表表5.3.1混凝土的抗混凝土的抗压压强强度度和养和养护时间护时间的的实测实测数据数据x23457912 14 17 21 28 56y35 42 47 53 59 65 68 73 76 82 86 995.3 数据拟合数据拟合5.4 插值和样条插值和样条 一、一、插插值Matlab提供了很多插提供了很多插值函数，函数，这些函数的使用方法些函数的使用方法大同小异，本大同小异，本节只介只介绍一一维插插值函数函数interp1，其他插，其他插值函数的用法函数的用法请读者参者参阅Matlab的关于的关于interp1的超文的超文

29、本帮助的本帮助的【See Also】条目。条目。曲线拟合是用带有噪声的曲线拟合是用带有噪声的曲线拟合是用带有噪声的曲线拟合是用带有噪声的“测量数据测量数据测量数据测量数据”构造出以某种方式最接近构造出以某种方式最接近构造出以某种方式最接近构造出以某种方式最接近“真实数据真实数据真实数据真实数据”的曲线方程，因为测量数据包含噪声，所以不要求拟合曲线的曲线方程，因为测量数据包含噪声，所以不要求拟合曲线的曲线方程，因为测量数据包含噪声，所以不要求拟合曲线的曲线方程，因为测量数据包含噪声，所以不要求拟合曲线穿越所有测量数据。穿越所有测量数据。穿越所有测量数据。穿越所有测量数据。插值和拟合不同，插值认为

30、原始数据是完全准确的，目的是用某种插值和拟合不同，插值认为原始数据是完全准确的，目的是用某种插值和拟合不同，插值认为原始数据是完全准确的，目的是用某种插值和拟合不同，插值认为原始数据是完全准确的，目的是用某种算法平滑的计算出这些原始数据之间的数据值。算法平滑的计算出这些原始数据之间的数据值。算法平滑的计算出这些原始数据之间的数据值。算法平滑的计算出这些原始数据之间的数据值。5.4 插值和样条插值和样条【调用格式用格式】yi=interp1(x,Y,xi,method)计算插算插值点自点自变量量为xi时的的值yi pp=interp1(x,Y,method,pp)用原始数据用原始数据获取插取插值

31、函数数据函数数据 yi=ppval(pp,xi)计算插算插值函数数据函数数据pp函数关系下的函数函数关系下的函数值【说明明】输入入变量量x，Y是原始数据向量是原始数据向量对。x的数据必的数据必须以以单调方式排列。方式排列。输入入变量量xi是插是插值点的自点的自变量坐量坐标向量。向量。输入入变量量method是插是插值方法，方法，Matlab可以可以选择的插的插值方法包括：方法包括：linear线性插性插值，缺省，缺省值。cublic三次多三次多项式插式插值 spline三次三次样条插条插值 nearst最最临近插近插值5.4 插值和样条插值和样条 输出出变量量yi为插插值点自点自变量量为xi时

32、的的计算算值。输出出变量量pp为插插值函数数据，里面保存着函数数据，里面保存着计算插算插值的的 表达式参数，用表达式参数，用 于描述相于描述相邻原始数据之原始数据之间的函数关系。的函数关系。可以通可以通过ppval函数函数计算算pp函数关系下的自函数关系下的自变量量xi的的 插插值结果果yi。例例5.4.1 已知已知1900年到年到1990年间，每隔年间，每隔10年美国的人口数量的统计数据年美国的人口数量的统计数据（单位：百万）依次为（单位：百万）依次为75.995，91.972，105.711，123.203，131.669，150.697，179.323，203.212，226.505，2

33、49.633，求在，求在1975年美国人年美国人口的数量，并绘制口的数量，并绘制1900到到1990年间每年的人口数量趋势图。年间每年的人口数量趋势图。5.4 插值和样条插值和样条二、样条二、样条 样条插值是常用的一种插值方法，其特点是精度高、最平滑，但是样条插值是常用的一种插值方法，其特点是精度高、最平滑，但是运算速度慢。经过样条插值后的曲线，除了在原始数据的端点外的其运算速度慢。经过样条插值后的曲线，除了在原始数据的端点外的其他数据点上都存在一阶和二阶导数，因此样条插值是非常平滑的。他数据点上都存在一阶和二阶导数，因此样条插值是非常平滑的。Matlab提供了专门用于样条的函数。提供了专门用

34、于样条的函数。【调用格式调用格式】yy=spline(x,y,xx)根据原始数据（根据原始数据（x，y）计算）计算xx的样条插值的样条插值yy pp=spline(x,y)根据原始数据（根据原始数据（x，y）计算分段）计算分段 样条函数数据样条函数数据pp dpp=fnder(pp)求求PP形式的样条函数的不定积分形式的样条函数的不定积分 ipp=fnint(pp)求求PP形式的样条函数的导数形式的样条函数的导数例例5.4.2 设函数设函数 ，用样条函数求用样条函数求和和只含有一个自变量的微分方程称为常微分方程（只含有一个自变量的微分方程称为常微分方程（ODE）。工程上的）。工程上的许多常微分

35、方程或者没有解析解，或者求解析解困难太大，这时可以选许多常微分方程或者没有解析解，或者求解析解困难太大，这时可以选择其数值解法。常微分方程分为初值问题和边值问题，本节只介绍初值择其数值解法。常微分方程分为初值问题和边值问题，本节只介绍初值问题的数值解法。问题的数值解法。5.5 常微分方程的数值解常微分方程的数值解一、一、ODE文件的编写格式文件的编写格式MATLAB中求解常微分方程的数值解是通过将其变形为一阶向量微中求解常微分方程的数值解是通过将其变形为一阶向量微分方程来实现的。用分方程来实现的。用MATLAB的的ODE解算指令解常微分方程，要编写表解算指令解常微分方程，要编写表示一阶向量微分

36、方程的函数示一阶向量微分方程的函数M文件，实现文件，实现 的微分的微分计算，其基本格式为：计算，其基本格式为：function DY=Fun(t,Y)其中：输入变量其中：输入变量t为时间变量，输入变量为时间变量，输入变量Y为列向量，为列向量，输出变量输出变量DY是是Y的一阶导数。的一阶导数。5.5 常微分方程的数值解常微分方程的数值解令令令令 则有则有则有则有 设常微分方程：设常微分方程：其初始条件为：其初始条件为：5.5 常微分方程的数值解常微分方程的数值解常微分方程化成一阶向量微分方程时，某些向量微分方程的向量解常微分方程化成一阶向量微分方程时，某些向量微分方程的向量解的各个分量的量级差别

37、较大，这对数值求解算法来说是很大的困难，这的各个分量的量级差别较大，这对数值求解算法来说是很大的困难，这种问题称之为种问题称之为刚性（刚性（stiff）问题）问题。MATLAB提供了很多常微分方程的解算函数，这些函数有些适用于提供了很多常微分方程的解算函数，这些函数有些适用于刚性方程，有些适用于非刚性方程，并且其使用的数值算法和解算精度刚性方程，有些适用于非刚性方程，并且其使用的数值算法和解算精度也各有不同，这些函数通称为也各有不同，这些函数通称为solve解算指令。解算指令。表表5.5.1中列出了各个解算指令的名称、精度和适用范围。中列出了各个解算指令的名称、精度和适用范围。二、二、solv

38、er解算指令解算指令5.5 常微分方程的数值解常微分方程的数值解solve指令指令解解题类题类型型精精 度度适用适用场场合合ode45非非刚刚性性一步法，一步法，4、5阶龙阶龙格格库库塔法，中等精度塔法，中等精度大多数大多数场场合的合的首首选选算法算法ode23非非刚刚性性一步法，一步法，2、3阶龙阶龙格格库库塔法，精度低塔法，精度低较较低精度低精度场场合，合，e-3ode113非非刚刚性性多步法，多步法，Adams算法，高低精度均可算法，高低精度均可ode45计计算算时间时间长时长时的替代的替代ode23t适度适度刚刚性性梯形法梯形法则则算法，精度低算法，精度低适度适度刚刚性性ode15s刚

39、刚性性多步法，中低精度多步法，中低精度ode45失失败时败时候适用候适用ode23s刚刚性性一步法，一步法，2阶阶Rosenbrock算式，精度低算式，精度低低精度低精度时时，比比ode15s有效有效ode23tb刚刚性性梯形法梯形法则则反向数反向数值值微分两微分两阶阶段法，段法，精度低精度低低精度低精度时时，比比ode15s有效有效表表5.5.1 solve解算指令解算指令5.5 常微分方程的数值解常微分方程的数值解【solve解算指令的解算指令的调用格式用格式】t,Y=solver(ODE_FUN,tspan,Y0,options)【说明说明】输入变量输入变量ODE_FUN是是ODE函数的

40、文件名。函数的文件名。输入变量输入变量tspan为求数值解的时间范围。当为求数值解的时间范围。当tspan=t0,tf时，表示求解时，表示求解 t0到到tf区间的数值解；当区间的数值解；当tspan=t0,t1,.,tf时，表示求解时，表示求解tspan指定指定 时间序列上的数值解。时间序列上的数值解。输入变量输入变量Y0是微分方程组的初始条件列向量。是微分方程组的初始条件列向量。输入变量输入变量options是可选择的算法参数，详见是可选择的算法参数，详见Matlab的帮助文件。的帮助文件。输出变量输出变量t是数值解的自变量数据列向量。是数值解的自变量数据列向量。输出变量输出变量Y是数值解。

41、是数值解。Y是一个矩阵，每一列代表了是一个矩阵，每一列代表了 向量解的一个分量在自变量向量解的一个分量在自变量t上的数值解。上的数值解。solver解算指令指的是表中的任意函数。解算指令指的是表中的任意函数。solver不仅可以求解常微分方程，还可以求解常微不仅可以求解常微分方程，还可以求解常微 分方程组。分方程组。5.5 常微分方程的数值解常微分方程的数值解例例5.5.1 求刚性常微分方程求刚性常微分方程 的解的解 其初始条件为其初始条件为 解：解：将微分方程的高阶导数写成低阶导数的组合，将微分方程的高阶导数写成低阶导数的组合，可以将原来的方程和初始条件转换为下面的一阶微分方程组和初始条件可

42、以将原来的方程和初始条件转换为下面的一阶微分方程组和初始条件例例5.5.2 求刚性常微分方程求刚性常微分方程 的解的解 其初始条件为其初始条件为 小小 结结 本章介绍常见的数值计算问题，包括线性代数、函数分析、数值本章介绍常见的数值计算问题，包括线性代数、函数分析、数值微积分、常微分方程等，重点介绍了如何利用微积分、常微分方程等，重点介绍了如何利用MATLAB提供的函数来提供的函数来实现数值计算，对于数学理论问题不做详细阐述，使用者很快就能应实现数值计算，对于数学理论问题不做详细阐述，使用者很快就能应用用MATLAB提供的强大函数进行复杂的方程组、微分、积分等运算。提供的强大函数进行复杂的方程

43、组、微分、积分等运算。MATLAB数值计算的结果为数值型数据，而不是数学上的解析表达式。数值计算的结果为数值型数据，而不是数学上的解析表达式。小小 结结51 用用LU分解求解下列分解求解下列线线性方程性方程组组。习习 题题52 设设 求函数求函数f在（在（0.5，0.5，0.5）附近的最小值。）附近的最小值。53 用一个多次三项式在区间用一个多次三项式在区间 内逼近函数内逼近函数 。习习 题题54 某某实验对实验对一根一根长长10m的的钢轨进钢轨进行行热热源的温度源的温度传传播播测试测试。用。用x表示表示测测量点距离（量点距离（m），用），用h表示表示测测量量时间时间（s），用），用T表示表示测测得各点的温得各点的温度（度（）,测测量量结结果如表所示。果如表所示。试用试用3次多项式插值求出在一分钟内钢管每隔次多项式插值求出在一分钟内钢管每隔10s、每隔、每隔0.5s处的温度。处的温度。55 设设有初有初值问题值问题：并与精确解作比并与精确解作比较较（精确解（精确解为为）。）。，试试求其数求其数值值解，解，