MATLAB编程与工程应用——第7章 MATLAB解方程与函数极值.ppt

1、MATLAB编程与工程应用,MATLAB解方程与函数极值,MATLAB编程与工程应用,教 师：曾晓东 电 话E_mail: QQ: 241747139,MATLAB解方程与函数极值,第7章 MATLAB解方程与函数极值,7.1 线性方程组求解 7.2 非线性方程数值求解 7.3 常微分方程初值问题的数值解法 7.4 函数极值,MATLAB解方程与函数极值,7.1 线性方程组求解,一、直接解法 对于线性方程组Ax=b，可以使用下面方法直接求解： 1利用左除运算符 x=Ab 2. 利用矩阵的逆 x=inv(A)*b 或 x=A(-1)*b 例7.1 用直接法求解线性方程组

2、 exp7_1.m,MATLAB解方程与函数极值,一、直接解法,3. 利用矩阵分解求解 矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesky分解、 Schur分解、奇异分解等。 对系数矩阵分解后再求解，可大大提高求解速度 (1) LU分解 矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。自变量A要求是方阵。 L,U=lu(A) x=U(Lb)=inv(L)*inv(U)*b A=L*U，其中L是一个下三角阵，U是一个上三角阵 L,U,P=lu(A) x=U(LP*b)=inv(L)*inv(U

3、)*P*b P*A=L*U，其中P是置换矩阵 例7.2 用LU分解求解例7.1 exp7_2.m,MATLAB解方程与函数极值,3. 利用矩阵分解求解,(2) QR分解 对矩阵X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。 QR分解只能对方阵进行。 Q,R=qr(A)x=R(Qb)=inv(Q)*inv(R)*b A=Q*R，其中Q是正交矩阵（即Q*Q=I)，R是上三角矩阵。 Q,R,E=qr(A)x=R(QE*b)=inv(R)*inv(Q)*E*b E*A=Q*R，其中E是置换矩阵。 例7.3 用QR分解求解例7.1 exp7_3.m,MATLAB解方程与函数极

4、值,7.1 线性方程组求解,二、迭代解法 迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。 1Jacobi迭代法 对于线性方程组Ax=b，如果A为非奇异方阵，即aii0， 则A可分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为： x=D-1(L+U)x+D-1b 与之对应的迭代公式为： x(k+1)=D-1(L+U)x(k)+D-1b 例7.4 利用jacobi迭代法求方程的解 exp7_4.m jacobi.m,MATLAB解方程

5、与函数极值,二、迭代解法,2Gauss-Serdel迭代法 在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到： x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些，稳定时间更快。 例7.5 利用Gauss-Serdel迭代法求方程的解 exp7_5.m gauseidel.m 例7.6 两种迭代法的比较 exp7_6.m,MATLAB解方程与函数极值,7.2 非线性方程数值求解,一

6、、单变量非线性方程求解 x=fzero(fname,x0,tol,trace) fname是待求根的函数文件名， x0为搜索的起点。 tol控制结果的相对精度，缺省时取tol=eps， trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省为0。 例7.7 求f(x)=x-10 x+2=0在x0=0.5附近的根。 funx.m exp7_7.m,MATLAB解方程与函数极值,7.2 非线性方程数值求解,二、非线性方程组求解 x=fsolve(fname,x0,option) fname用于定义需求解的非线性方程组的函数文件名。 x0为搜索的起点。 option优化选项设定 最优

7、化工具箱提供了20多个选项，用户可以使用optimget命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。 Display选项决定函数调用时中间结果的显示方式，其中off为不显示，iter表示每步都显示，final只显示最终结果。optimset(Display,off)将设定Display选项为off。,MATLAB解方程与函数极值,7.2 非线性方程数值求解,例7.8 求下列非线性方程组在(0.5,0.5) 附近的数值解。,(1)建立函数文件func.m表示非线性方程组 (2)调用fsolve求方程组的解exp7_8.m (3)回代以验证结果func(x

8、),MATLAB解方程与函数极值,7.3 常微分方程初值问题的数值解法,一、龙格库塔法简介 考虑求如下常微分方程（组）初值问题的数值解,常用方法有： 1. Euler(欧拉)方法：ym+1=ym+hf(xm,ym) 这是一个递推方法，从x0处开始递推求解，h为步长,MATLAB解方程与函数极值,一、龙格库塔法简介,2. Runge-Kutta(龙格-库塔）法 递推公式的形式为：,这是一种线性多步法，其中h为步长，s为阶次，bi为待定系数,如s=2时,中点公式,MATLAB解方程与函数极值,一、龙格库塔法简介,如s=3时,如s=4时,Kutta三阶算法,古典龙格-库塔法,MATLAB解方程与函数

9、极值,一、龙格库塔法简介,3. 刚性问题 在用微分方程描述的一个变化过程中，若往往又包含着多个相互作用但变化速度相差十分悬殊的子过程，这样一类过程就认为具有“刚性”。 描述这类过程的微分方程初值问题称为“刚性问题”。 刚性问题广泛存在于化学反应过程、电力系统、航空航天、热核反应等领域 例如，宇航飞行器自动控制系统一般包含两个相互作用但效应速度相差十分悬殊的子系统，一个是控制飞行器质心运动的系统，当飞行器速度较大时，质心运动惯性较大，因而相对来说变化缓慢；另一个是控制飞行器运动姿态的系统，由于惯性小，相对来说变化很快，因而整个系统就是一个刚性系统。,MATLAB解方程与函数极值,7.3 常微分方

10、程初值问题的数值解法,二、MATLAB用于求常微分方程数值解的函数 ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb,MATLAB解方程与函数极值,二、MATLAB用于求常微分方程数值解的函数,一般调用格式（以ode45为例） t,y=ode45(fname,tspan,y0) 其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。 tspan形式为t0,tf,表示求解区间。要获得问题在时间点t0,t1,tf上的解，则令tspan=t0,t1,tf（要求是单调的） y0是初始状态列向量。 t和y分别给出时间向量和相应的状态向量。 例7.9 求解描述振荡器的经典的Ver der Pol微分方程 verderpol.m exp7_9.m,MATLAB解方程与函数极值,7.4 函数极值,1. 求单变量函数的最小值点 x=fminbnd(fname,x1,x2) 2. 求多变量函数的最小值点 x=f