matlab总复习课件.ppt

1、科学计算与科学计算与MATLABMATLABn1 1 科学计算科学计算数值计算的定义数值计算的定义数值计算的研究过程数值计算的研究过程数值计算的误差分析数值计算的误差分析数值计算的稳定性与收敛性数值计算的稳定性与收敛性n2 MATLAB2 MATLAB简介简介MATLABMATLAB的类别的类别MATLABMATLAB的概述的概述n3 3 课程的要求与学习方法课程的要求与学习方法课程的内容课程的内容课程的要求课程的要求课程的学习与考核课程的学习与考核n教材与主要参考书教材与主要参考书n小结小结n如果存在一适当小的如果存在一适当小的正数正数r ，使得，使得则称则称r为为相相对误差限对误差限。 例

2、：例：x=15, (x) =2; r(x)=2/15=13.33%y=1000,(y)=5 , r(y)=5/1000=0.5%rrxxxxxexe)()(2.有效数字：有效数字：定义：若近似值定义：若近似值x的误差限是某一位的半个的误差限是某一位的半个单位，该位到单位，该位到x的第一位非零数字共有的第一位非零数字共有n位，位，就说就说x有有n位有效数字，位有效数字，x可表示为可表示为其中，其中，a10.且且 a1, a2, , an 都是都是 0 9 中的中的任一整数。任一整数。其绝对误差限满足：其绝对误差限满足：mnaaax10. 021nmxxxe1021)(例例. 测量一物体的长度为测

3、量一物体的长度为954cm，问测量数据的相，问测量数据的相对误差限多大？对误差限多大？解解:因实际问题所截取的近似数，其绝对误差限一因实际问题所截取的近似数，其绝对误差限一般不超过最小刻度的半个单位，故当般不超过最小刻度的半个单位，故当x=954cm时，时，有有(x)0.5cm，而而x的相对误差的相对误差er(x) 0.5/954=0.0005241.0.00053=0.053 %故故r(x)=0.053 %.第二讲第二讲 基础知识（数基础知识（数值、符号计算）值、符号计算）n1 1 MATLABMATLAB的启动与运行的启动与运行MatlabMatlab的启动和退出的启动和退出MatlabM

4、atlab的工作窗口的工作窗口MatlabMatlab的指令行的操作的指令行的操作MatlabMatlab的帮助系统的帮助系统n2 MATLAB2 MATLAB的矩阵与数值计算的矩阵与数值计算MatlabMatlab数据类型数据类型MatlabMatlab的变量及其命名的变量及其命名MatlabMatlab矩阵及其运算矩阵及其运算MatlabMatlab的数组关系的数组关系/ /逻辑运算逻辑运算MatlabMatlab的多项式运算的多项式运算n3 MATLAB3 MATLAB的符号计算的符号计算n小结小结1.5 Matlab的帮助系统的帮助系统MATLAB帮助命令包括帮助命令包括help、lo

5、okfor以及模糊查询以及模糊查询1 help命令命令 在命令窗口中输入在命令窗口中输入help命令将会显示当前帮助系统中所命令将会显示当前帮助系统中所包含的所有项目，即搜索路径中所有的目录名称。也可包含的所有项目，即搜索路径中所有的目录名称。也可以通过以通过help加函数名来显示该函数的帮助说明。加函数名来显示该函数的帮助说明。2 lookfor命令命令 对搜索范围内的对搜索范围内的M文件的第一行进行关键字搜索，条件文件的第一行进行关键字搜索，条件比较宽松。若对比较宽松。若对M文件进行全文搜索，加上文件进行全文搜索，加上-all选项。选项。3 模糊查询模糊查询 MATLAB 6.0以上的版本

6、提供了一种类似模糊查询的以上的版本提供了一种类似模糊查询的命令查询方法，只需要输入命令的前几个字母，然后按命令查询方法，只需要输入命令的前几个字母，然后按Tab键，系统就会列出所有以这几个字母开头的命令。键，系统就会列出所有以这几个字母开头的命令。数值显示格式数值显示格式 任何任何MATLAB的语句的执行结果都可以在屏幕的语句的执行结果都可以在屏幕上显示，同时赋值给指定的变量，没有指定变上显示，同时赋值给指定的变量，没有指定变量时，赋值给一个特殊的变量量时，赋值给一个特殊的变量ans，数据的显，数据的显示格式由示格式由format命令命令控制。控制。 format只是影响结果的显示，不影响其计

7、算与只是影响结果的显示，不影响其计算与存储；存储；MATLAB总是以双字长浮点数（双精总是以双字长浮点数（双精度）来执行所有的运算。度）来执行所有的运算。 如果结果为整数，则显示没有小数；如果结果不如果结果为整数，则显示没有小数；如果结果不是整数，则输出形式有：是整数，则输出形式有：format (short)：短格式（：短格式（5位定点数）位定点数）99.1253format long：长格式（：长格式（15位定点数位定点数 ） 99.12345678900000format short e：短格式：短格式e方式方式 9.9123e+001format long e：长格式：长格式e方式方式

8、 9.912345678900000e+001format bank：2位十进制位十进制 99.12format hex：十六进制格式：十六进制格式符号变量和符号表达式的创建符号变量和符号表达式的创建 MATLAB提供了两个建立符号量的命令：提供了两个建立符号量的命令：syms 和和 sym 说明：说明：1. a1 a2 a3.为需要定义为符号量的标识符，不能是数字、函数表为需要定义为符号量的标识符，不能是数字、函数表达式或方程式。达式或方程式。2. a1 a2 a3.均不用引号界定，而用均不用引号界定，而用空格空格分开。分开。1 用用syms创建符号量创建符号量 格式：格式：syms a1

9、a2 a3说明：说明：1. 输入参量输入参量A可以是数字、字符串也可以是字符串变量名、字符表可以是数字、字符串也可以是字符串变量名、字符表达式或字符方程。达式或字符方程。2. 回车后将回车后将A定义成符号量并赋值给标识符。定义成符号量并赋值给标识符。2 用用sym创建符号量、符号表达式创建符号量、符号表达式 格式：标识符格式：标识符=sym（A）3.3.矩阵算法矩阵算法 和和 数组算法数组算法矩阵算法矩阵算法 把矩阵看作一个整体，各种运算完全按照线性代数代表把矩阵看作一个整体，各种运算完全按照线性代数代表的矩阵运算法则进行，运算的书写形式和运算符号都与的矩阵运算法则进行，运算的书写形式和运算符

10、号都与矩阵理论完全相同。矩阵理论完全相同。数组算法数组算法把矩阵看作由其元素构成的一组数据（数组），各种运把矩阵看作由其元素构成的一组数据（数组），各种运算是在参与运算矩阵的对应元素之间进行的数与数的运算是在参与运算矩阵的对应元素之间进行的数与数的运算，这种运算方便对大批数据的处理和一次求出多个函算，这种运算方便对大批数据的处理和一次求出多个函数值。数值。 数组算法的运算符主要有数组算法的运算符主要有 . .* * ./ . . ./ . . Matlab运算符运算符 算术运算符算术运算符操作符操作符说说 明明+加加 A+B A B必须大小相同，或一个是标量必须大小相同，或一个是标量-减减 A

11、-B A B必须大小相同，或一个是标量必须大小相同，或一个是标量*矩阵乘矩阵乘 A*B A 的列数等于的列数等于B的行数，一个可以是标量的行数，一个可以是标量.*数组乘数组乘 A.*B A B必须大小相同，一个可以是标量必须大小相同，一个可以是标量矩阵左除矩阵左除 AB =A-1* B 等效于等效于A*X=B求求X inv(A).数组左除数组左除 A.B Bij/Aij/矩阵右除矩阵右除 A/B =A*B-1 等效于等效于X*B=A求求X./数组右除数组右除 A./B Aij/Bij矩阵乘方矩阵乘方 Ap A自乘自乘p次次.数组乘方数组乘方 A.p A中每个元素的中每个元素的p次方次方如：如：

12、a=1 2;3 4；b= 3 5; 5 9c=a+b d=a-bc= d= 4 7 -2 -3 8 13 -2 -5a*b=13 23; 29 51a/b=-0.50 0.50;3.50 1.50ab=-1 -1;2 3a3=37 54; 81 118a.*b=3 10;15 36a./b=0.33 0.40;0.60 0.44a.b=3.00 2.50;1.67 2.25a.3= 1 8; 27 64 只有维数相同的矩阵才能进只有维数相同的矩阵才能进行加减运算。行加减运算。 注意只有当两个矩阵中前一注意只有当两个矩阵中前一个矩阵的列数和后一个矩阵个矩阵的列数和后一个矩阵的行数相同时，才可以进

13、行的行数相同时，才可以进行乘法运算。乘法运算。ab运算等效于运算等效于求求a*x=b的解；而的解；而a/b等效于等效于求求x*b=a的解。只有方阵才的解。只有方阵才可以求幂。可以求幂。 转置：对于实矩阵用（转置：对于实矩阵用（）符号或（）符号或（.）求转置结果）求转置结果是一样的；然而对于含复数的矩阵，则（是一样的；然而对于含复数的矩阵，则（）将同）将同时对复数进行共轭处理，而时对复数进行共轭处理，而 （.）则只是将其排列形）则只是将其排列形式进行转置。式进行转置。4.常用矩阵运算常用矩阵运算逆矩阵与行列式计算逆矩阵与行列式计算求逆求逆：inv(A)；求行列式求行列式：det(A)要求矩阵必须

14、为方阵要求矩阵必须为方阵5、矩阵分解、矩阵分解（1）奇异值分解）奇异值分解U,S,V=svd(A)例：例：a = 9 8 6 8可以验证：可以验证：u*u=Iv*v=Iu*s*v=a求矩阵求矩阵A的奇异值及分解矩阵，满的奇异值及分解矩阵，满足足U*S*V=A，其中，其中U、V矩阵为矩阵为正交矩阵（正交矩阵（U*U=I），），S矩阵为对矩阵为对角矩阵，它的对角元素即角矩阵，它的对角元素即A矩阵的矩阵的奇异值。奇异值。u,s,v=svd(a)u = 0.7705 -0.6375 0.6375 0.7705s = 15.5765 0 0 1.5408v = 0.6907 -0.7231 0.7231

15、 0.6907（2）特征值分解）特征值分解V,D=eig(A)例：例： a = 9 8 6 8v,d=eig(a)v = 0.7787 -0.7320 0.6274 0.6813d = 15.4462 0 0 1.5538求矩阵求矩阵A的特征向量的特征向量V及特征值及特征值D，满足满足A*V=V*D。其中。其中D的对角线的对角线元素为特征值，元素为特征值，V的列为对应的特的列为对应的特征向量。如果征向量。如果D=eig(A)则只返回则只返回特征值。特征值。（3）正交分解）正交分解Q,R=qr(A)例：例： a = 9 8 6 8q,r=qr(a)q = -0.8321 -0.5547 -0.5

16、547 0.8321r = -10.8167 -11.0940 0 2.2188将矩阵将矩阵A做正交化分解，使得做正交化分解，使得Q*R=A，其中，其中Q为正交矩阵（其为正交矩阵（其范数为范数为1，指令，指令norm(Q)=1)，R为为对角化的上三角矩阵。对角化的上三角矩阵。（4）三角分解）三角分解L,U=lu(A)将将A做对角线分解，使得做对角线分解，使得A=L*U,其中其中L为下三角矩阵，为下三角矩阵，U为上三角矩阵。为上三角矩阵。注意：注意：L实际上是一个实际上是一个“心理上心理上”的下三角矩阵，它事实的下三角矩阵，它事实上是一个置换矩阵上是一个置换矩阵P的逆矩阵与一个真正下三角矩阵的逆

17、矩阵与一个真正下三角矩阵L1（其对角线元素为（其对角线元素为1）的乘积。）的乘积。L1,U1,P=lu(A)例：例：a=1 2 3;4 5 6;7 8 9 比较：比较： l1,u1,p=lu(a) l,u=lu(a)3 3、MATLABMATLAB的符号计算的符号计算 A = sym(A) sum(A) sum(A) det(A) inv(A) eig(A) svd(A)符号计算常规：符号计算常规：1 下面的MATLAB语句中正确的有： a) 2api; b) record_1=3+4i c) a=2.0, d) c=1+6j 2. 已知水的黏度随温度的变化公式如下，其中a=0.03368，b

18、=0.000221，计算温度t为20，30，40度时的粘度分别是： 021atbt0为0水的黏度， 值为31.785 10； a、 b为常数， 分别为0.03368、 0.000221。 （a）0.0018 0.0010 0.0007 (b) 0.0010 0.0007 0.0005 (c) 1.7850e-003 1.0131e-003 6.6092e-004 (d) 1.0131e-003 6.6092e-004 4.6772e-004 第三讲第三讲 基础知识（图基础知识（图形、程序设计）形、程序设计）n1 1 MATLABMATLAB的图形功能的图形功能二维图形的基本函数二维图形的基本函

19、数三维图形的基本函数三维图形的基本函数图形窗口函数图形窗口函数简单动画函数简单动画函数n2 MATLAB2 MATLAB的程序设计的程序设计M M文件介绍文件介绍MatlabMatlab的程序结构的程序结构MatlabMatlab的文件操作的文件操作MatlabMatlab的函数操作的函数操作n小结小结1.1.1 基本二维图形基本二维图形函数函数1：plot：绘制二维数据图形：绘制二维数据图形格式：格式：1. plot(X,S)2. plot(X,Y,S)3. plot(X1,Y1,S1,X2,Y2,S2 ,.,X3,Y3,S3)说明说明:参数参数S控制数据点的标记、曲线类型和曲线色彩，三者控

20、制数据点的标记、曲线类型和曲线色彩，三者置于一对单引号内。置于一对单引号内。常用的绘图选项常用的绘图选项选选 项项含含 义义选选 项项含含 义义-实线实线*用星号标出数据点用星号标出数据点-虚线虚线.用点号标出数据点用点号标出数据点 :点线点线o用圆圈号标出数据点用圆圈号标出数据点-.点划线点划线x用叉号标出数据点用叉号标出数据点b蓝色蓝色+用加号标出数据点用加号标出数据点g绿色绿色s用正方形标出数据点用正方形标出数据点r红色红色D用菱形标出数据点用菱形标出数据点c青色青色V用下三角标出数据点用下三角标出数据点m洋红洋红 用上三角标出数据点用上三角标出数据点y黄色黄色用右三角标出数据点用右三角

21、标出数据点w白色白色H用六角形标出数据点用六角形标出数据点P用五角形标出数据点用五角形标出数据点函数函数2： fplot：绘制：绘制y=f(x)图形图形格式：格式：fplot(fname,lims, S)说明说明:其中其中fname为函数名或单引号界定的函数表达式，为函数名或单引号界定的函数表达式，lims为为x,y的取值范围，的取值范围，S定义与定义与plot函数相同。函数相同。函数函数3： ezplot：绘制隐函数图形：绘制隐函数图形格式：格式：ezplot(f, xmin,xmax,ymin,ymax)说明说明:在区间在区间xminxxmax和和yminyymax绘制绘制f(x,y) =

22、 0的图的图形。形。,默认区间默认区间-2x2和和-2y21 f(x,y)=0格式：格式：ezplot(X,Y, tmin,tmax)说明说明:在区间在区间tmin t tmax绘制绘制x=X(t)和和y=Y(t)的图形的图形,默认区默认区间间0t22 x=X(t) y=Y(t)格式：格式：ezplot(f, a,b)说明说明:在区间在区间axb绘制绘制y= f(x)的图形的图形,默认区间默认区间-2x x=1:5; y=1.5 1.8 4 3.4 5.7; a,b=LZXEC(x,y)a = 1.0000b = 0.2800 常用的线性变换常用的线性变换 函数函数变换后的函数变换后的函数)(

23、xfy BAxykxaey akxylnlnbxaybxay1cbxayacxaby12)(baxcybaxcy21baxcykxakxbycln)ln(kxaxey akxxylnln4、matlab拟合函数拟合函数线性拟合函数线性拟合函数格式：格式：p=linefit(x,y)说明：说明：x，y 输入同维数据向量输入同维数据向量p 输出拟合多项式的系数向量输出拟合多项式的系数向量多项式拟合函数多项式拟合函数格式：格式：p=polyfit(x,y,m)说明：说明：x，y 输入同维数据向量输入同维数据向量m 拟合多项式次数拟合多项式次数p 输出拟合多项式的系数输出拟合多项式的系数 向量向量4.

24、 设y=span1,x,x2，用最小二乘法拟合如下表数据。x0.51.01.52.02.53.0y1.75 2.45 3.81 4.80 8.00 8.60计算的结果为：a) 0.4900 1.2501 0.8560b) 0.8560 1.2501 0.4900c) -0.6341 3.8189 -3.7749d) 3.8189 -3.7749 2.8533 x=0.5:0.5:3.0; y=1.75,2.45,3.81,4.80,8.00,8.60; a=polyfit(x,y,2)第六讲数值积分与微分第六讲数值积分与微分n引言引言n数值积分数值积分矩形积分近似计算矩形积分近似计算梯形积分近

25、似计算梯形积分近似计算抛物线形积分近似计算抛物线形积分近似计算牛顿牛顿-科茨（科茨（Newton-Cotes）公式）公式自适应（自适应（Simpson）求积法）求积法高斯高斯(Gauss)求积法求积法n数值微分数值微分nMATLAB的积分和微分函数的积分和微分函数n小结小结MATLAB中提供的积分函数：中提供的积分函数：矩形积分近似计算矩形积分近似计算梯形积分近似计算梯形积分近似计算抛物线形积分近似计算抛物线形积分近似计算牛顿牛顿-科茨（科茨（Newton-Cotes）公式）公式高斯高斯(Gauss)求积法求积法自适应（自适应（Simpson）求积法）求积法龙贝格积分法龙贝格积分法样条函数求积

26、分法样条函数求积分法简单奇异积分简单奇异积分重积分计算重积分计算 4.2 matlab数值计算定积分数值计算定积分格式：格式：Z = trapz(x,y)（1） 梯形求积公式梯形求积公式说明说明x,y分别代表数目相同分别代表数目相同 的阵列或矩阵。的阵列或矩阵。y与与x具有一定的函数关系。具有一定的函数关系。格式：格式：I,n=quad(fname,a,b,tol,trace)（2）自适应辛卜生法）自适应辛卜生法说明说明返回参数返回参数I即定积分值，即定积分值，n为被积函数的调用次数。为被积函数的调用次数。fname是被积函数名。定义被积函数须用数值运算符。是被积函数名。定义被积函数须用数值运

27、算符。a和和b分别是定积分的下限和上限。分别是定积分的下限和上限。tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非控制是否展现积分过程，若取非0则展现积分过程，则展现积分过程，取取0则不展现，缺省时取则不展现，缺省时取trace=0。格式：格式：I,n=quad8(fname,a,b,tol,trace) （1） Newton-Cotes法法说明说明tol用来控制积分精度，缺省时取用来控制积分精度，缺省时取tol=10-6。其它与。其它与quad相同。相同。该函数可以更精确地求出定积分的值，且计算精度相该函数可以更精确地求出定

28、积分的值，且计算精度相同条件下，函数调用的步数明显小于同条件下，函数调用的步数明显小于quad函数，从而函数，从而保证能以更高的效率求出所需的定积分值。保证能以更高的效率求出所需的定积分值。格式：格式：I,n=quadl(fname,a,b,tol,trace) （4）高斯）高斯-洛巴托（洛巴托（ Gauss-Labatto）法）法格式：格式：I=dblquad(fname,xmin,xmax,ymin,ymax,tol ,trace) （5） 二重积分二重积分 4.3 matlab数值计算微分数值计算微分格式：格式：f =diff(fun)F=diff(fun,x)F=diff(fun,x,

29、n)F=diff(S)F=diff(S,x)F=diff(S,x,n)（1） 数值微分和差分数值微分和差分一元函数求导一元函数求导格式：格式：gx,gy=gradient(F)gx,gy=gradient(F,H)nx,ny,nz=surfnorm(X,Y,Z)二元、多元函数求导二元、多元函数求导clear all;df1 = diff(sin(x);df2 = diff(sqrt(x)df3 = diff(log(x)df4= diff(log(x*y),x)df5 = diff(sin(x*y),y)df6 = diff(sin(x*y),x,3)Ans:df1 = cos(x) df2

30、= 1/2/x(1/2) df3 = 1/x df4 = 1/x df5 = cos(x*y)*x df6 = -cos(x*y)*y3clear all; F= 1 1.2 1.4 2.3 5;0 -0.6 3 4 2;-1 7 7.2 9 1.4;gx,gy=gradient(F) %计算梯度计算梯度n=surfnorm(F) %计算法向量计算法向量6. 用Newton-Cotes方法计算如下积分 12532xdxx （a）133.6625 (b) 23.8600 (c) 87.9027 (d) -1.6180 function f=fun(x)f=x.*x.*sqrt(2*x+3)qua

31、dl(fun,1,3,1e-10) 或quadl(exp(-x/2),1,3,1e-10)或fun=(x)exp(-x/2); quadl(fun,1,3,1e-10)syms xf = ( x 3 + x 2 + x + 1 ) ( 1 / 3 ) -sqrt(x2+x+1)*log(exp(x)+x)/xlimit(f,x,inf)7. y=ln(1+x)，求x=1时y的近似值。 a) -0.25 b) 0.5 c) -0.6137 d) -1.6137 syms xy=log(1+x)f=diff(y,2)subs(f,1)第七讲常微分方程数值解法第七讲常微分方程数值解法n引言引言n欧拉

32、近似方法欧拉近似方法n龙格龙格-库塔库塔(R-K)方法方法nMATLAB的常微分方程函数的常微分方程函数n小结小结4、 MATLAB的常微分方程函数的常微分方程函数ode45 ode23 ode113 ode15s ode23s ode23t ode23tb格式格式x,y=ode45( fun, x0,xn, y0,option说明：说明：适用于求解一阶常微分方程组适用于求解一阶常微分方程组fun定义微分方程组的函数文件名定义微分方程组的函数文件名x0,xn求解区域求解区域 y0初始条件向量初始条件向量option可选参数，由可选参数，由ODESET函数设置，比较复杂函数设置，比较复杂x输出自

33、变量向量，输出自变量向量，y输出输出y, y , y ,.,.2 , 1 ),()(iyxfxyiii没有一种算法可以有效地解决所有的没有一种算法可以有效地解决所有的 ODE 问题，因此问题，因此MATLAB 提供了多种提供了多种ODE函数。函数。函数函数ODE类类型型特点特点说明说明ode45非刚性非刚性单步法；单步法；4，5 阶阶 R-K 方法；方法；累计截断误差为累计截断误差为 (x)3大部分场合的首选方法大部分场合的首选方法ode23非刚性非刚性单步法；单步法；2，3 阶阶 R-K 方法；方法；累计截断误差为累计截断误差为 (x)3使用于精度较低的情形使用于精度较低的情形ode113非

34、刚性非刚性多步法；多步法；Adams算法；高低精算法；高低精度均可到度均可到 10-310-6计算时间比计算时间比 ode45 短短ode23t适度刚性适度刚性 采用梯形算法采用梯形算法适度刚性情形适度刚性情形ode15s刚性刚性多步法；多步法；Gears 反向数值微分；反向数值微分；精度中等精度中等若若 ode45 失效时，可尝失效时，可尝试使用试使用ode23s刚性刚性单步法；单步法；2 阶阶Rosebrock 算法；算法；低精度低精度当精度较低时，计算时当精度较低时，计算时间比间比 ode15s 短短ode23tb刚性刚性梯形算法；低精度梯形算法；低精度当精度较低时，计算时当精度较低时，

35、计算时间比间比ode15s短短例题：用例题：用MATLAB的符号解法的符号解法,求解常微分方程：求解常微分方程：23xdyxyxedx dsolve(Dy+3*x*y=x*exp(-x2) ans =(1/3*exp(-x*(x-3*t)+C1)*exp(-3*x*t) dsolve(Dy+3*x*y=x*exp(-x2),x) ans = exp(-x2)+exp(-3/2*x2)*C1例题：采用例题：采用ODE45求解描述某刚性问题的微分方程：求解描述某刚性问题的微分方程：123121323123(0)0,(0)10.51(0)1yy yyyy yyyy yy function dy =

36、rigid(t,y)dy = zeros(3,1); % 行向量行向量dy(1) = y(2) * y(3);dy(2) = -y(1) * y(3);dy(3) = -0.51 * y(1) * y(2);options = odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-5);T,Y = ode45(rigid,0 12,0 1 1,options);plot(T,Y(:,1),-,T,Y(:,2),-.,T,Y(:,3),.)legend(y1,y2,y3)例题：用例题：用MATLAB函数函数ode23,ode45,ode113,求解多阶常微分方程求解多阶常

37、微分方程：323232233(1)1,(1)10, (1)30,1,10 xd yd ydyxedxdxdxyyyx12232332333233xdyydxdyydxdyeyydxxxx1223333010000103230 xydYdyYdxdxyexxxfunction dy=myfun03(x,y)dy=zeros(3,1) %初始化变量初始化变量dydy(1)=y(2); %dy(1)表示表示y的一阶导数的一阶导数,其等于其等于y的第二列值的第二列值dy(2)=y(3); %dy(2)表示表示y的二阶导数的二阶导数dy(3)=2*y(3)/x3+3*y(2)/x3+3*exp(2*x)

38、/x3 %dy(3)表示表示y的三阶导数的三阶导数% 用用ode23 ode45 ode113解多阶微分方程解多阶微分方程clear,clcx23,y23=ode23(myfun03,1,10,1 10 30);x45,y45=ode45(myfun03,1,10,1 10 30);x113,y113=ode113(myfun03,1,10,1 10 30);figure(1) %第一幅图第一幅图plot(x23,y23(:,1),*r,x45,y45(:,1),ob,x113,y113(:,1),+g) %作出各种函数所得结果作出各种函数所得结果legend(ode23解解,ode45解解,

39、ode113解解)title(ODE函数求解结果函数求解结果)figure(2)plot(x45,y45) %以以ode45为例作出函数以及其各阶导数图为例作出函数以及其各阶导数图legend(y,y一阶导数一阶导数,y两阶导数两阶导数)title(y,y一阶导数一阶导数,y二阶导数函数图二阶导数函数图)10. 采用ODE45求解如下多阶常微分方程，并求出当x1.8505时的函数值。 a) 31.6441 b) 74.6907 c) 118.7862 d) 63.2564 323232233(1)1,(1)10, (1)30,1,10 xd yd ydyxedxdxdxyyyx建立求解函数文件

40、myfun03function dy=myfun03(x,y)dy=zeros(3,1) %初始化变量初始化变量dy，该行可以没有，该行可以没有dy(1)=y(2); %dy(1)表示表示y的一阶导数的一阶导数,其等于其等于y的第二列值的第二列值dy(2)=y(3); %dy(2)表示表示y的二阶导数的二阶导数dy(3)=2*y(3)/x3+3*y(2)/x3+3*exp(2*x)/x3 %dy(3)表示表示y的三阶导数的三阶导数求解过程：x45,y45=ode45(myfun03,1,10,1 10 30);第八讲非线性方程求根第八讲非线性方程求根n引言引言n二分法二分法n迭代法迭代法nNe

41、wton迭代法迭代法nMATLAB的非线性方程求根函数的非线性方程求根函数n小结小结Matlab中中zeroin的算法实现是的算法实现是fzero. x = fzero(FUN,x0) %x0可以是数，或区间可以是数，或区间x = fzero(FUN,x0,options)x,fval= fzero(FUN,x0,options)x,fval,exitflag = fzero(FUN,x0,options)5、MATLAB的非线性方程求根函数的非线性方程求根函数例：例： 求方程在求方程在x=0.5附近的根。附近的根。210 xx x=fzero(x2+x-1,0.5)x = 0.6180例：采

42、用牛顿下山法求方程在区间例：采用牛顿下山法求方程在区间1.2，2上的一个上的一个根。根。220 xx x=fzero(sqrt(x)-x3+2,0 2)x = 1.4759 r=NewtonDown(sqrt(x)-x3+2,1.2,2)r = 1.47595. 求方程在x=0.5附近的根. 21xx a) 0.6180 b) -1.1719e-25 c) 1 d) -1.6180 第九讲解线性方程组的直接解法第九讲解线性方程组的直接解法n引言引言nGauss消元法消元法n列主元素消元法列主元素消元法n误差分析误差分析nMATLAB的线性方程组求解函数的线性方程组求解函数1n小结小结2、Gau

43、ss消元法消元法n2.1 基本思想：基本思想：逐步消去未知元，将方程组化为与逐步消去未知元，将方程组化为与其等价的上三角方程组求解。其等价的上三角方程组求解。n2.2 分两步：分两步： 第一步第一步: 消元过程消元过程，将方程组消元化为等价的上三角形，将方程组消元化为等价的上三角形方程组方程组; 第二步第二步: 回代过程回代过程， 解上三角形方程组解上三角形方程组,得原方程组的得原方程组的解。解。 )1()1()1(2)1(22)1(2211212111nnnnnnnnnnbxabxaxabxaxaxa2.3 Gauss消元的目的消元的目的a11x1+ a12x2+ a1nxn = b1a21

44、x1+ a22x2+ a2nxn = b2an1x1+ an2x2+ annxn = bn原始方程组原始方程组约化方程组约化方程组2.4.1 消元过程消元过程(化一般方程组为上三角方程组化一般方程组为上三角方程组)以四阶为例：以四阶为例：其系数增广矩阵为：其系数增广矩阵为：4444343242141343433323213124243232221211414313212111bxaxaxaxabxaxaxaxabxaxaxaxabxaxaxaxa444434241334333231224232221114131211baaaabaaaabaaaabaaaaA第一轮消元：第一轮消元：计算计算3个

45、数个数: m21 m31 m41T = a21 a31 a41T / a11 用用-m21乘矩阵第一行后加到矩阵第二行乘矩阵第一行后加到矩阵第二行; 用用-m31乘矩阵第一行后加到矩阵第三行乘矩阵第一行后加到矩阵第三行; 用用-m41乘矩阵第一行后加到矩阵第四行乘矩阵第一行后加到矩阵第四行;其系数增广矩阵变为：其系数增广矩阵变为：)()()()()()()()()()()()()(1414414314213134133132121241231221141312111baaabaaabaaabaaaaA)0(11a第二轮消元：第二轮消元：计算计算2个数：个数：m32 m42T = a32(1)

46、a42(1)T / a22(1) 用用-m32乘矩阵第二行后加到矩阵第三行乘矩阵第二行后加到矩阵第三行; 用用-m42乘矩阵第二行后加到矩阵第四行乘矩阵第二行后加到矩阵第四行;其系数增广矩阵变为：其系数增广矩阵变为：)()()()()()()()()()()(2424424323234233121241231221141312112baabaabaaabaaaaA)0()1(22a第三轮消元：第三轮消元：计算计算: m43=a43(2)/a33(2)用用-m43乘矩阵第三行后加到矩阵第四行乘矩阵第三行后加到矩阵第四行;其系数增广矩阵变为：其系数增广矩阵变为：)()()()()()()()()(

47、)(3434423234233121241231221141312113babaabaaabaaaaA)0()2(33a其对应的上三角方程组为其对应的上三角方程组为)()()()()()()()()(3443442342343233124124312321221414313212111bxabxaxabxaxaxabxaxaxaxa 若对于一般的线性方程组若对于一般的线性方程组Ax=b，其消元过程的计算，其消元过程的计算公式为公式为： （k=1,2,n-1）), 1(), 1,(), 1()1()1()()1()1()()1(nkibmbbnkjiamaankiaamkkikkikikkjik

48、kijkijkkkikik）（记ijijaa)0(2.4.2 回代过程回代过程(解上三角方程组解上三角方程组)上三角方程组的一般形式为：上三角方程组的一般形式为： 其中其中a11ann0nnnnnnnnbxabxaxabxaxaxa2222211212111) 1 , 2 , 1()(1niaxabxabxiinijjijiinnnn回代过程的计算公式回代过程的计算公式：定理：定理： 约化的主元素约化的主元素ak+1,k+1(k) 0 (k=0,1,n-1)的充分必要条件是的充分必要条件是 矩阵矩阵A的各阶顺序主子式的各阶顺序主子式不为零。即不为零。即01111 kkkkkaaaaD注：注：对

49、角线上的元素对角线上的元素ak+1,k+1(k)在在Gauss消去法中作消去法中作用突出，称约化的主元素。用突出，称约化的主元素。推论推论: 如果如果A的顺序主子式的顺序主子式Dk 0 (k=1,n-1)，则则Gauss消元法中的约化主元可以表示为消元法中的约化主元可以表示为) 1,.,3 , 2(/1)(1, 1111nkDDaDakkkkk例例 用高斯消元法求解方程组用高斯消元法求解方程组 32563034253432321xxx 32303452536432202230445 . 0064324200445 . 006432 446245 . 0432321xxxx1= -13, x2

50、= 8, x3 = 2m21=3/2m31=4/2m32=-3/0.5 clear all;A=2 3 4;3 5 2;4 3 30;b=6 5 32;L,U=lu(A);x=U(Lb)x = -13.0000 8.0000 2.0000MATLAB 处理：处理：线性代数方程组线性代数方程组11 11221331121 1222233221 12233nnnnnnnnnnnA xA xA xA xBA xA xA xA xBA xA xA xA xB可以用矩阵形式表示为可以用矩阵形式表示为111213111212223222123 nnnnnnnnnAAAAxBAAAAxBxBAAAA即：即：