第7章+MATLAB解方程與函數(shù)極值（MATLAB程序設(shè)計(jì)教程電子教案）.ppt

第7章MATLAB解方程與函數(shù)極值7 1線性方程組求解7 2非線性方程數(shù)值求解7 3常微分方程初值問題的數(shù)值解法7 4函數(shù)極值 7 1線性方程組求解7 1 1直接解法1 利用左除運(yùn)算符的直接解法對(duì)于線性方程組Ax b 可以利用左除運(yùn)算符 求解 x A b 例7 1用直接解法求解下列線性方程組 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 x A b 2 利用矩陣的分解求解線性方程組矩陣分解是指根據(jù)一定的原理用某種算法將一個(gè)矩陣分解成若干個(gè)矩陣的乘積 常見的矩陣分解有LU分解 QR分解 Cholesky分解 以及Schur分解 Hessenberg分解 奇異分解等 1 LU分解矩陣的LU分解就是將一個(gè)矩陣表示為一個(gè)交換下三角矩陣和一個(gè)上三角矩陣的乘積形式 線性代數(shù)中已經(jīng)證明 只要方陣A是非奇異的 LU分解總是可以進(jìn)行的 MATLAB提供的lu函數(shù)用于對(duì)矩陣進(jìn)行LU分解 其調(diào)用格式為 L U lu X 產(chǎn)生一個(gè)上三角陣U和一個(gè)變換形式的下三角陣L 行交換 使之滿足X LU 注意 這里的矩陣X必須是方陣 L U P lu X 產(chǎn)生一個(gè)上三角陣U和一個(gè)下三角陣L以及一個(gè)置換矩陣P 使之滿足PX LU 當(dāng)然矩陣X同樣必須是方陣 實(shí)現(xiàn)LU分解后 線性方程組Ax b的解x U L b 或x U L Pb 這樣可以大大提高運(yùn)算速度 例7 2用LU分解求解例7 1中的線性方程組 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 L U lu A x U L b 或采用LU分解的第2種格式 命令如下 L U P lu A x U L P b 2 QR分解對(duì)矩陣X進(jìn)行QR分解 就是把X分解為一個(gè)正交矩陣Q和一個(gè)上三角矩陣R的乘積形式 QR分解只能對(duì)方陣進(jìn)行 MATLAB的函數(shù)qr可用于對(duì)矩陣進(jìn)行QR分解 其調(diào)用格式為 Q R qr X 產(chǎn)生一個(gè)一個(gè)正交矩陣Q和一個(gè)上三角矩陣R 使之滿足X QR Q R E qr X 產(chǎn)生一個(gè)一個(gè)正交矩陣Q 一個(gè)上三角矩陣R以及一個(gè)置換矩陣E 使之滿足XE QR 實(shí)現(xiàn)QR分解后 線性方程組Ax b的解x R Q b 或x E R Q b 例7 3用QR分解求解例7 1中的線性方程組 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 Q R qr A x R Q b 或采用QR分解的第2種格式 命令如下 Q R E qr A x E R Q b 3 Cholesky分解如果矩陣X是對(duì)稱正定的 則Cholesky分解將矩陣X分解成一個(gè)下三角矩陣和上三角矩陣的乘積 設(shè)上三角矩陣為R 則下三角矩陣為其轉(zhuǎn)置 即X R R MATLAB函數(shù)chol X 用于對(duì)矩陣X進(jìn)行Cholesky分解 其調(diào)用格式為 R chol X 產(chǎn)生一個(gè)上三角陣R 使R R X 若X為非對(duì)稱正定 則輸出一個(gè)出錯(cuò)信息 R p chol X 這個(gè)命令格式將不輸出出錯(cuò)信息 當(dāng)X為對(duì)稱正定的 則p 0 R與上述格式得到的結(jié)果相同 否則p為一個(gè)正整數(shù) 如果X為滿秩矩陣 則R為一個(gè)階數(shù)為q p 1的上三角陣 且滿足R R X 1 q 1 q 實(shí)現(xiàn)Cholesky分解后 線性方程組Ax b變成R Rx b 所以x R R b 例7 4用Cholesky分解求解例7 1中的線性方程組 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 R chol A Errorusing cholMatrixmustbepositivedefinite命令執(zhí)行時(shí) 出現(xiàn)錯(cuò)誤信息 說明A為非正定矩陣 7 1 2迭代解法迭代解法非常適合求解大型系數(shù)矩陣的方程組 在數(shù)值分析中 迭代解法主要包括Jacobi迭代法 Gauss Serdel迭代法 超松弛迭代法和兩步迭代法 1 Jacobi迭代法對(duì)于線性方程組Ax b 如果A為非奇異方陣 即aii 0 i 1 2 n 則可將A分解為A D L U 其中D為對(duì)角陣 其元素為A的對(duì)角元素 L與U為A的下三角陣和上三角陣 于是Ax b化為 x D 1 L U x D 1b與之對(duì)應(yīng)的迭代公式為 x k 1 D 1 L U x k D 1b這就是Jacobi迭代公式 如果序列 x k 1 收斂于x 則x必是方程Ax b的解 Jacobi迭代法的MATLAB函數(shù)文件Jacobi m如下 function y n jacobi A b x0 eps ifnargin 3eps 1 0e 6 elseifnargin epsx0 y y B x0 f n n 1 end 例7 5用Jacobi迭代法求解下列線性方程組 設(shè)迭代初值為0 迭代精度為10 6 在命令中調(diào)用函數(shù)文件Jacobi m 命令如下 A 10 1 0 1 10 2 0 2 10 b 9 7 6 x n jacobi A b 0 0 0 1 0e 6 2 Gauss Serdel迭代法在Jacobi迭代過程中 計(jì)算時(shí) 已經(jīng)得到 不必再用 即原來的迭代公式Dx k 1 L U x k b可以改進(jìn)為Dx k 1 Lx k 1 Ux k b 于是得到 x k 1 D L 1Ux k D L 1b該式即為Gauss Serdel迭代公式 和Jacobi迭代相比 Gauss Serdel迭代用新分量代替舊分量 精度會(huì)高些 Gauss Serdel迭代法的MATLAB函數(shù)文件gauseidel m如下 function y n gauseidel A b x0 eps ifnargin 3eps 1 0e 6 elseifnargin epsx0 y y G x0 f n n 1 end 例7 6用Gauss Serdel迭代法求解下列線性方程組 設(shè)迭代初值為0 迭代精度為10 6 在命令中調(diào)用函數(shù)文件gauseidel m 命令如下 A 10 1 0 1 10 2 0 2 10 b 9 7 6 x n gauseidel A b 0 0 0 1 0e 6 例7 7分別用Jacobi迭代和Gauss Serdel迭代法求解下列線性方程組 看是否收斂 命令如下 a 1 2 2 1 1 1 2 2 1 b 9 7 6 x n jacobi a b 0 0 0 x n gauseidel a b 0 0 0 7 2非線性方程數(shù)值求解7 2 1單變量非線性方程求解在MATLAB中提供了一個(gè)fzero函數(shù) 可以用來求單變量非線性方程的根 該函數(shù)的調(diào)用格式為 z fzero fname x0 tol trace 其中fname是待求根的函數(shù)文件名 x0為搜索的起點(diǎn) 一個(gè)函數(shù)可能有多個(gè)根 但fzero函數(shù)只給出離x0最近的那個(gè)根 tol控制結(jié)果的相對(duì)精度 缺省時(shí)取tol eps trace 指定迭代信息是否在運(yùn)算中顯示 為1時(shí)顯示 為0時(shí)不顯示 缺省時(shí)取trace 0 例7 8求f x x 10 x 2 0在x0 0 5附近的根 步驟如下 1 建立函數(shù)文件funx m functionfx funx x fx x 10 x 2 2 調(diào)用fzero函數(shù)求根 z fzero funx 0 5 z 0 3758 7 2 2非線性方程組的求解對(duì)于非線性方程組F X 0 用fsolve函數(shù)求其數(shù)值解 fsolve函數(shù)的調(diào)用格式為 X fsolve fun X0 option 其中X為返回的解 fun是用于定義需求解的非線性方程組的函數(shù)文件名 X0是求根過程的初值 option為最優(yōu)化工具箱的選項(xiàng)設(shè)定 最優(yōu)化工具箱提供了20多個(gè)選項(xiàng) 用戶可以使用optimset命令將它們顯示出來 如果想改變其中某個(gè)選項(xiàng) 則可以調(diào)用optimset 函數(shù)來完成 例如 Display選項(xiàng)決定函數(shù)調(diào)用時(shí)中間結(jié)果的顯示方式 其中 off 為不顯示 iter 表示每步都顯示 final 只顯示最終結(jié)果 optimset Display off 將設(shè)定Display選項(xiàng)為 off 例7 9求下列非線性方程組在 0 5 0 5 附近的數(shù)值解 1 建立函數(shù)文件myfun m functionq myfun p x p 1 y p 2 q 1 x 0 6 sin x 0 3 cos y q 2 y 0 6 cos x 0 3 sin y 2 在給定的初值x0 0 5 y0 0 5下 調(diào)用fsolve函數(shù)求方程的根 x fsolve myfun 0 5 0 5 optimset Display off x 0 63540 3734 將求得的解代回原方程 可以檢驗(yàn)結(jié)果是否正確 命令如下 q myfun x q 1 0e 009 0 23750 2957可見得到了較高精度的結(jié)果 7 3常微分方程初值問題的數(shù)值解法7 3 1龍格 庫塔法簡(jiǎn)介7 3 2龍格 庫塔法的實(shí)現(xiàn)基于龍格 庫塔法 MATLAB提供了求常微分方程數(shù)值解的函數(shù) 一般調(diào)用格式為 t y ode23 fname tspan y0 t y ode45 fname tspan y0 其中fname是定義f t y 的函數(shù)文件名 該函數(shù)文件必須返回一個(gè)列向量 tspan形式為 t0 tf 表示求解區(qū)間 y0是初始狀態(tài)列向量 t和y分別給出時(shí)間向量和相應(yīng)的狀態(tài)向量 例7 10設(shè)有初值問題 試求其數(shù)值解 并與精確解相比較 精確解為y t 1 建立函數(shù)文件funt m functionyp funt t y yp y 2 t 2 4 t 1 2 求解微分方程 t0 0 tf 10 y0 2 t y ode23 funt t0 tf y0 求數(shù)值解y1 sqrt t 1 1 求精確解t y y1 y為數(shù)值解 y1為精確值 顯然兩者近似 例7 11求解著名的VanderPol方程 例7 12有Lorenz模型的狀態(tài)方程 試?yán)L制系統(tǒng)相平面圖 7 4函數(shù)極值MATLAB提供了基于單純形算法求解函數(shù)極值的函數(shù)fmin和fmins 它們分別用于單變量函數(shù)和多變量函數(shù)的最小值 其調(diào)用格式為 x fmin fname x1 x2 x fmins fname x0 這兩個(gè)函數(shù)的調(diào)用格式相似 其中fmin函數(shù)用于求單變量函數(shù)的最小值點(diǎn) fname是被最小化的目標(biāo)函數(shù)名 x1和x2限定自變量的取值范圍 fmins函數(shù)用于求多變量函數(shù)的最小值點(diǎn) x0是求解的初