数值计算 精选课件.ppt

《数值计算 精选课件.ppt》由会员分享，可在线阅读，更多相关《数值计算 精选课件.ppt（44页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

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

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

4、e-14数量数量级）例例例例5.1.25.1.2第五页，本课件共有44页三、奇异三、奇异值分解分解1.矩矩阵的奇异的奇异值分解分解任意矩任意矩阵A可可进行奇异行奇异值分解，即存在酉矩分解，即存在酉矩阵U和和V，使下面等式成立，使下面等式成立其中其中 称称为矩矩阵的奇异的奇异值。5.1 线性代数线性代数【调用格式用格式】s=svd(A)求矩求矩阵A的奇异的奇异值，s为由奇异由奇异值构成的向量构成的向量U,S,V=svd(A)矩矩阵A的奇异的奇异值分解分解第六页，本课件共有44页5.1 线性代数线性代数2.矩矩阵结构特征的奇异构特征的奇异值描述描述 矩矩阵的奇异的奇异值可以描述矩可以描述矩阵的的结

5、构特征。有关矩构特征。有关矩阵结构特征的构特征的MATLAB 函数有如下几种。函数有如下几种。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)在指定容差在指定容差t

6、ol下，求矩下，求矩阵A的的广广义逆逆，tol可以省略。可以省略。第七页，本课件共有44页5.1 线性代数线性代数 四、四、线性方程性方程组的解的解形如形如 的的线性方程性方程组中中独立方程的个数独立方程的个数等于等于独立未知参数的个数称独立未知参数的个数称为恰定方程；恰定方程；独立方程的个数独立方程的个数大于大于独立未知参数的个数称独立未知参数的个数称为超定方程；超定方程；独立方程的个数独立方程的个数小于小于独立未知参数的个数称独立未知参数的个数称为欠定方程。欠定方程。第八页，本课件共有44页5.1 线性代数线性代数1.左除运算符法左除运算符法 对于于一般的非奇异矩一般的非奇异矩阵A，可以求

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

8、件共有44页5.2 函数分析函数分析 一、函数的零点一、函数的零点1.多多项式的根式的根 通通过roots函数来求取多函数来求取多项式全部的根。式全部的根。【调用格式用格式】r=roots(p)多多项式求根函数，求取多式求根函数，求取多项式的全部根式的全部根【说明明】p为多多项式的式的系数行向量系数行向量，r为多多项式所有式所有根构成的列向量根构成的列向量第十一页，本课件共有44页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区区间内的零点。内的零点。第十二页，本课件共有44页5.2 函数分析函数分

10、析输入入变量量options为优化迭代化迭代选项，是一个，是一个结构体。构体。输出出变量量x为零点零点处的自的自变量量值，输出出变量量fval为零点零点处的函数的函数值。输出出变量量exitflag表示函数中止表示函数中止计算的条件。算的条件。若若exitflag0表示找到零点后退出。表示找到零点后退出。例例5.2.1 求函数求函数在在区区间内的所有零点。内的所有零点。第十三页，本课件共有44页5.2 函数分析函数分析【调用格式用格式】x,fval,exitflag,output=fsolve(fun,x0,options)完整格式完整格式x=fsolve(fun,x0)多元非多元非线性方程求

11、解的最性方程求解的最简格式格式【说明明】fun表示自表示自变量量为向量的函数，可以是字符串、内向量的函数，可以是字符串、内联函数、函数句柄。函数、函数句柄。输入入变量量x0为零点的初始猜零点的初始猜测向量（自向量（自变量量值）。）。输出出变量量x为零点零点处的自的自变量量值，输出出变量量fval为零点零点处的函数的函数值。例例5.2.2 求方程求方程组在在附近的解。附近的解。3.多元函数的零点多元函数的零点多元函数的精确零点可以用多元函数的精确零点可以用fsolve函数求取，必函数求取，必须提供零点的大致位置才能提供零点的大致位置才能进行数行数值搜索。搜索。第十四页，本课件共有44页5.2 函

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

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

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

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

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

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

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表示表示选用的用的积分算法，可以省略，确是算法的分算法，可以省略，确是算法的 是是quad，还可以可以选则quadl或者用或者用户

19、自定自定义的的积分分 算法算法M函数文件的函数句柄。函数文件的函数句柄。本函数无法本函数无法计算内重算内重积分上下限分上下限为函数表达式的情况。函数表达式的情况。5.2 函数分析函数分析第二十页，本课件共有44页求形如求形如 的的闭型三重数型三重数值积分的函数分的函数为triplequad。【调用格式用格式】q=triplequad(fun,x1,x2,y1,y2,z1,z2,tol,method)3.三重数三重数值积分分【说明明】fun为被被积函数函数 ，可以用字符串、内，可以用字符串、内联函数和函数函数和函数M文件的文件的函数句函数句 柄表示。柄表示。积分分变量用量用x，y，z表示，表示，

20、从内到外从内到外积分分变量依次量依次为x，y，z。x2和和x1是是变量量x的的积分上下限，分上下限，y2和和y1是是变量量y的的积分分 上下限，上下限，z2和和z1是是变量量z的的积分上下限。分上下限。其他事其他事项同函数同函数dblquad。5.2 函数分析函数分析例例5.2.6 计算定算定积分分第二十一页，本课件共有44页5.3 数据拟合数据拟合已知某些离散的原始数据，建立一个曲已知某些离散的原始数据，建立一个曲线方程，方程，让它以最佳的方它以最佳的方式反映原始数据的式反映原始数据的变化化趋势，尽量避免出，尽量避免出现局部的波局部的波动，这种方法称种方法称为拟合。不要求合。不要求拟合曲合曲

21、线经过所有的原始数据。最佳方式通常是指用所有的原始数据。最佳方式通常是指用拟合得到的数学模型合得到的数学模型计算出来的算出来的计算数据和原始数据之算数据和原始数据之间的的误差的平差的平方和最小，方和最小，这种方式称种方式称为最小二乘。通常是最小二乘。通常是已知曲已知曲线方程方程的形式（数的形式（数学模型的学模型的结构），根据原始数据来构），根据原始数据来计算曲算曲线方程的参数方程的参数（数学模型的（数学模型的参数）。本参数）。本节介介绍Matlab中如何中如何实现最小二乘最小二乘拟合。合。第二十二页，本课件共有44页5.3 数据拟合数据拟合一、多一、多项式式拟合合多多项式式拟合是用一个合是用一

22、个n阶多多项式模型去式模型去拟合原始数据的方法，需要合原始数据的方法，需要计算的算的模型参数包括多模型参数包括多项式的式的阶次和多次和多项式的系数。式的系数。Matlab提供了提供了polyfit函数函数实现多多项式式拟合。合。【调用格式用格式】p=polyfit(x,y,n)根据根据给定的数据定的数据(x,y)，计算算n阶拟合多合多项式的式的系数向量系数向量p ye=polyval(p,x)计算自算自变量量为x时多多项式式p的的值（估（估计值）第二十三页，本课件共有44页 5.3 数据拟合数据拟合 【说明明】多多项式式拟合一般不要超合一般不要超过5阶，否，否则计算算误差差变大。大。拟合多合多

23、项式只在原始数据范式只在原始数据范围内可以保内可以保证精度，超出范精度，超出范围使用使用 拟合多合多项式无法保式无法保证预报的精度。的精度。例例5.3.1：已知五个数据点：已知五个数据点：1,5.5,2,43,3,128,4,290,5,498，试画出画出这五个点五个点拟合的三次曲合的三次曲线。第二十四页，本课件共有44页二、最小二乘二、最小二乘拟合合1.线性最小二乘性最小二乘拟合合对于于线性数学模型的参数估性数学模型的参数估计，可以用形如，可以用形如Y=Ax+b的一的一阶多多项式式拟合来估合来估计参数。某些非参数。某些非线性模型性模型经过变量替量替换也可以也可以转换为线性模型，也性模型，也可

24、以采用可以采用线性估性估计方法。方法。【调用格式用格式】a=lsqlin(Fun,a0,lb,ub,options,p1,p2,.)a=lsqnonlin(Fun,a0)5.3 数据拟合数据拟合第二十五页，本课件共有44页例例5.3.2：对于非于非线性数学模型性数学模型 ，其中，其中a和和b为模型模型参数，参数，采用采用变量替量替换将其将其转换为线性数学模型。性数学模型。解：解：对已知数学模型取自然已知数学模型取自然对数有数有 令令 则有的有的线性数学模型性数学模型5.3 数据拟合数据拟合第二十六页，本课件共有44页函数的功能是求取向量的函数的功能是求取向量的值，使向量函数，使向量函数 满足足

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

26、页，本课件共有44页输入入变量量a0为向量，表示所求向量，表示所求变量量a的初始猜的初始猜测值。输入入变量量lb表示所求表示所求变量量a的下限，的下限，输入入变量量ub表示所求表示所求变量量a的上的上限。如果没有上下限限制，可以用空矩限。如果没有上下限限制，可以用空矩阵表示，如果没有下限表示，如果没有下限可以用可以用-inf表示。表示。从第六个从第六个输入入变量开始的参数量开始的参数p1,p2等是向等是向Fun函数函数传递的参数，的参数，其名字其名字应该和定和定义Fun函数的函数的输入入变量名相一致。量名相一致。输出出变量量a表示所求表示所求变量的量的值。如果将如果将 构造构造为残差函数，且残

27、差函数，且 为被估被估计参数参数向量，向量，则lsqnonlin就就变为求解非求解非线性最小二乘估性最小二乘估计的函数。的函数。5.3 数据拟合数据拟合第二十八页，本课件共有44页例例5.3.3 混凝土的抗混凝土的抗压强度度 随养随养护时间 的延的延长而增加，其数学模型而增加，其数学模型为现将一批混凝土作成将一批混凝土作成12个个测试块，记录了养了养护时间x（日）及抗（日）及抗压强度度y （kg/cm2）的数据，如表）的数据，如表5.3.1 所示。所示。试估估计该数学模型的参数数学模型的参数a1，a2，a3，a4。表表5.3.1混凝土的抗混凝土的抗压强度度和养和养护时间的的实测数据数据x234

28、57912 14 17 21 28 56y35 42 47 53 59 65 68 73 76 82 86 995.3 数据拟合数据拟合第二十九页，本课件共有44页5.4 插值和样条插值和样条 曲曲曲曲线线拟拟合合合合是用是用是用是用带带有噪声的有噪声的有噪声的有噪声的“测测量数据量数据量数据量数据”构造出以某种方式最构造出以某种方式最构造出以某种方式最构造出以某种方式最接近接近接近接近“真真真真实实数据数据数据数据”的曲的曲的曲的曲线线方程，因方程，因方程，因方程，因为测为测量数据包含噪声，所以量数据包含噪声，所以量数据包含噪声，所以量数据包含噪声，所以不不不不要求要求要求要求拟拟合曲合曲合

29、曲合曲线线穿越所有穿越所有穿越所有穿越所有测测量数据量数据量数据量数据。插插插插值值和和和和拟拟合不同，合不同，合不同，合不同，插插插插值认为值认为原始数据是完全准确的原始数据是完全准确的原始数据是完全准确的原始数据是完全准确的，目的是，目的是，目的是，目的是用某种算法平滑的用某种算法平滑的用某种算法平滑的用某种算法平滑的计计算出算出算出算出这这些原始数据之些原始数据之些原始数据之些原始数据之间间的数据的数据的数据的数据值值。第三十页，本课件共有44页5.4 插值和样条插值和样条【调用格式用格式】yi=interp1(x,Y,xi,method)计算插算插值点自点自变量量为xi时的的值yi p

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

31、值 nearst最最临近插近插值一、一、插插值第三十一页，本课件共有44页5.4 插值和样条插值和样条 例例5.4.1 已知已知1900年到年到1990年年间，每隔，每隔10年美国的人口数量的年美国的人口数量的统计数据（数据（单位：百万）依次位：百万）依次为75.995，91.972，105.711，123.203，131.669，150.697，179.323，203.212，226.505，249.633，求在，求在1975年美国人口的数量，并年美国人口的数量，并绘制制1900到到1990年年间每年的人口数量每年的人口数量趋势图。t=1900:10:1990;p=75.995 91.972

32、 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633;interp1(t,p,1975)interp1(t,p,2000)%无法无法计算自算自变量超出范量超出范围的插的插值，返回，返回值为NaNpp=interp1(t,p,spline,pp);%求取人口的插求取人口的插值函数数据函数数据tt=1900:1990;yy=ppval(pp,tt);%根据插根据插值函数数据函数数据计算插算插值结果果plot(t,p,o,tt,yy,b)%绘制插制插值函数的曲函数的曲线第三十二页，本课件共有44页5.4 插值和样条插值和样条

33、二、二、样条条 样条插条插值是常用的一种插是常用的一种插值方法，其特点是方法，其特点是精度高、最平滑精度高、最平滑，但是运算，但是运算速速度慢度慢。经过样条插条插值后的曲后的曲线，除了在原始数据的端点外的其他数据点上都，除了在原始数据的端点外的其他数据点上都存存在一在一阶和二和二阶导数数，因此，因此样条插条插值是非常平滑的。是非常平滑的。Matlab提供了提供了专门用于用于样条条的函数。的函数。【调用格式用格式】yy=spline(x,y,xx)根据原始数据（根据原始数据（x，y）计算算xx的的样条插条插值yy pp=spline(x,y)根据原始数据（根据原始数据（x，y）计算分段算分段 样

34、条函数数据条函数数据pp dpp=fnder(pp)求求PP形式的形式的样条函数的不定条函数的不定积分分 ipp=fnint(pp)求求PP形式的形式的样条函数的条函数的导数数第三十三页，本课件共有44页例例5.4.2 设函数函数 ，用用样条函数求条函数求和和5.4 插值和样条插值和样条dx=0.05;%设定采样时间间隔设定采样时间间隔x=(0:dx:1)*pi;y=sin(x).*cos(x);%求采样点数据求采样点数据pp=spline(x,y);%求求f(x)的样条函数的样条函数ipp=fnint(pp);%求求pp样条函数的不定积分，近似为样条函数的不定积分，近似为f(x)的不定积分的

35、不定积分dpp=fnder(pp);%求求pp样条函数的微分，近似为样条函数的微分，近似为f(x)的微分的微分ppval(ipp,1,2)*-1;1%ppval(ipp,2)-ppval(ipp,1)，求定积分，求定积分ppval(dpp,2)第三十四页，本课件共有44页只含有一个自只含有一个自变量的微分方程称量的微分方程称为常微分方程（常微分方程（ODE）。工程上的）。工程上的许多多常微分方程或者没有解析解，或者求解析解困常微分方程或者没有解析解，或者求解析解困难太大，太大，这时可以可以选择其其数数值解法解法。5.5 常微分方程的数值解常微分方程的数值解一、一、ODE文件的文件的编写格式写格

36、式MATLAB中求解常微分方程的数中求解常微分方程的数值解是通解是通过将其将其变形形为一一阶向量微分方程向量微分方程来来实现的。要的。要编写表示写表示一一阶向量微分方程向量微分方程的函数的函数M文件，文件，实现 的微分的微分计算，其基本格式算，其基本格式为：function DY=Fun(t,Y)其中：其中：输入入变量量t为时间变量，量，输入入变量量Y为列向量，列向量，输出出变量量DY是是Y的一的一阶导数。数。第三十五页，本课件共有44页5.5 常微分方程的数值解常微分方程的数值解令令令令 设常微分方程：常微分方程：其初始条件其初始条件为：设常微分方程：常微分方程：其初始条件其初始条件为：令令

37、第三十六页，本课件共有44页5.5 常微分方程的数值解常微分方程的数值解常微分方程化成一常微分方程化成一阶向量微分方程向量微分方程时，某些向量微分方程的向量解，某些向量微分方程的向量解的各个分量的量的各个分量的量级差差别较大，大，这对数数值求解算法来求解算法来说是很大的困是很大的困难，这种种问题称之称之为刚性（性（stiff）问题。MATLAB提供了很多常微分方程的解算函数，提供了很多常微分方程的解算函数，这些函数有些适用于些函数有些适用于刚性方性方程，有些适用于非程，有些适用于非刚性方程，并且其使用的数性方程，并且其使用的数值算法和解算精度也各有不同，算法和解算精度也各有不同，这些函数通称些

38、函数通称为solve解算指令。解算指令。表表5.5.1中列出了各个解算指令的名称、精度和适用范中列出了各个解算指令的名称、精度和适用范围。二、二、solver解算指令解算指令第三十七页，本课件共有44页5.5 常微分方程的数值解常微分方程的数值解solve指令解题类型精 度适用场合ode45非刚性一步法，4、5阶龙格库塔法，中等精度大多数场合的首选算法ode23非刚性一步法，2、3阶龙格库塔法，精度低较低精度场合，e-3ode113非刚性多步法，Adams算法，高低精度均可ode45计算时间长时的替代ode23t适度刚性梯形法则算法，精度低适度刚性ode15s刚性多步法，中低精度ode45失败

39、时候适用ode23s刚性一步法，2阶Rosenbrock算式，精度低低精度时，比ode15s有效ode23tb刚性梯形法则反向数值微分两阶段法，精度低低精度时，比ode15s有效表表5.5.1 solve解算指令解算指令第三十八页，本课件共有44页5.5 常微分方程的数值解常微分方程的数值解【solve解算指令的解算指令的调用格式用格式】t,Y=solver(ODE_FUN,tspan,Y0,options)【说明明】输入入变量量ODE_FUN是是ODE函数的文件名。函数的文件名。输入入变量量tspan为求数求数值解的解的时间范范围。当tspan=t0,tf时，表示求解 0到tf区间的数值解；

40、当tspan=t0,t1,.,tf时，表示求解tspan指定 时间序列上的数值解。输入入变量量Y0是微分方程是微分方程组的初始条件列向量。的初始条件列向量。输入入变量量options是可是可选择的算法参数，的算法参数，详见Matlab的帮助文件。的帮助文件。输出出变量量t是数是数值解的自解的自变量数据列向量。量数据列向量。输出出变量量Y是数是数值解。解。Y是一个矩阵，每一列代表了向量解的一个分量在自变量t上的数值解。solver解算指令指的是表中的任意函数。解算指令指的是表中的任意函数。solver不不仅可以求解常微分方程，可以求解常微分方程，还可以求解常微分方程可以求解常微分方程组。第三十九

41、页，本课件共有44页5.5 常微分方程的数值解常微分方程的数值解例例5.5.1 求求刚性常微分方程性常微分方程 的解的解 其初始条件其初始条件为 解：解：将微分方程的高将微分方程的高阶导数写成低数写成低阶导数的数的组合，合，可以将原来的方程和初始条件可以将原来的方程和初始条件转换为下面的一下面的一阶微分方程微分方程组和初始条件和初始条件第四十页，本课件共有44页5.5 常微分方程的数值解常微分方程的数值解例例5.5.2 求求刚性常微分方程性常微分方程 的解的解 其初始条件其初始条件为 第四十一页，本课件共有44页 小小结 本章介本章介绍常常见的数的数值计算算问题，包括，包括线性代数、函数分析、

42、数性代数、函数分析、数值微微积分、常微分方程等，重点介分、常微分方程等，重点介绍了如何利用了如何利用MATLAB提供的函数来提供的函数来实现数数值计算，算，对于数学理于数学理论问题不做不做详细阐述，使用者很快就能述，使用者很快就能应用用MATLAB提供的提供的强大函数大函数进行复行复杂的方程的方程组、微分、微分、积分等运算。分等运算。MATLAB数数值计算的算的结果果为数数值型数据，而不是数学上的解析表达式。型数据，而不是数学上的解析表达式。小小 结结第四十二页，本课件共有44页51 用用LU分解求解下列分解求解下列线性方程性方程组。习习 题题52 设 求函数求函数f在（在（0.5，0.5，0.5）附近的最小）附近的最小值。53 用一个多次三用一个多次三项式在区式在区间 内逼近函数内逼近函数 。第四十三页，本课件共有44页感谢大家观看第四十四页，本课件共有44页