代數(shù)方程和常微分方程求解.ppt

第8章代數(shù)方程和常微分方程求解 代數(shù)方程是未知數(shù)和常數(shù)進(jìn)行有限次代數(shù)運(yùn)算所組成的方程 它包括有理方程和無理方程 代數(shù)方程的解稱為的根或零點(diǎn) 其求解一般是通過代數(shù)幾何來進(jìn)行 微分方程是含有一個或是多個導(dǎo)數(shù)的方程 只有一個自變量及其導(dǎo)數(shù)的微分方程稱為常微分方程 包含兩個以上自變量及其偏導(dǎo)數(shù)的微分方程稱為偏微分方程 工程上許多物理規(guī)律 設(shè)計(jì)過程的模擬和評價(jià) 凡是涉及質(zhì)量和能量運(yùn)動設(shè)計(jì)分析的問題 都最終歸結(jié)到微分方程 8 1代數(shù)方程求解 8 1 1代數(shù)方程圖解法符號繪圖函數(shù)fplot 和ezplot 也可以用于圖解法求代數(shù)方程的根 它適用于求解維數(shù)較少的一維方程或二維方程組 對于一維方程圖解 其解就是函數(shù)曲線與x軸交點(diǎn)所對應(yīng)的變量數(shù)值 如果有多個交點(diǎn) 則表示該方程有多個解 如果沒有交點(diǎn) 則表示該方程沒有解 例如 在例5 3使用符號繪圖函數(shù)繪制代數(shù)方程的圖形 圖5 3左圖 中可見 函數(shù)在區(qū)間 5 5 內(nèi)與x軸有3個交點(diǎn) 因此該代數(shù)方程該區(qū)間內(nèi)有3個實(shí)根 對于二維方程組圖解 其解就是兩條函數(shù)曲線的交點(diǎn)所對應(yīng)的坐標(biāo)數(shù)值 如果只有1個交點(diǎn) 或切點(diǎn) 則表示該方程組有1個解 如果有2個交點(diǎn) 則表示該方程組有2個解 如果沒有交點(diǎn) 則表示該方程沒有解 例8 1用圖解法求解二維聯(lián)立方程 a 2 b 2 定義橫軸區(qū)間ezplot x 2 y 2 1 69 a b axis equal 控制坐標(biāo)軸比例相等holdon gridon ezplot 2 4 x 3 y 1 5 a b line a b 0 0 line 0 0 b a xlabel bfx ylabel bfy title bf二維代數(shù)方程組的圖解法 gtext bff 1 x 2 y 2 1 3 2 gtext bff 2 2 4x 3 y 1 5 8 1 2代數(shù)方程的解析解求非線性方程或方程組解析解的函數(shù)調(diào)用格式 X solve fun x 其中 fun是符號方程的函數(shù)表達(dá)式 x是自變量 X是解析解 應(yīng)當(dāng)指出 函數(shù)solve fun x 也可以用于求線性方組的解析解 例8 2求非線性解方程組解析解 二維非線性方程組的解析解symsabxy f1 x 2 x y a f2 y 2 x y b disp 二維非線性方程組的解析解 X Y solve f1 f2 x y M文件運(yùn)行結(jié)果 二維非線性方程組的解析解 x a a b 1 2 a a b 1 2 Y 1 a b 1 2 b 1 a b 1 2 b 8 1 3線性方程組的數(shù)值解最簡便方法是使用矩陣左除或是矩陣求逆的方法 求解線性方程組AX b X A bX inv A b其中 A是方程組的系數(shù)矩陣 b是常數(shù)向量 X是解析解 例8 3求線性方程組的數(shù)值解 線性方程組的數(shù)值解AA 1 1 1 1 1 2 1 4 2 3 1 5 3 1 2 11 bb 5 2 2 0 線性方程組常數(shù)向量disp 采用矩陣左除求出線性方程組的解 xx AA bbdisp 采用矩陣求逆求出線性方程組的解 zx inv AA bbdisp 計(jì)算殘量 r AA zx bbdisp 計(jì)算殘量的模 R norm r M文件運(yùn)行結(jié)果 采用矩陣左除或矩陣求逆求出線性方程組的解 xx zx 1 00002 00003 0000 1 0000計(jì)算殘量 r 1 0e 014 0 08880 2220 0 44410 1776計(jì)算殘量的模 R 5 3475e 015 8 1 4非線性方程的數(shù)值解1 一維非線性方程對于一維非線性方程求解 可以看作是單變量的極小化問題 通過不斷縮小搜索區(qū)間來逼近一維問題的真解 因此 可以使用一維非線性方程優(yōu)化解函數(shù)來求解 其調(diào)用格式是 x fx flag fzero fun x0 其中 輸入?yún)?shù)中 fun是非線性方程的函數(shù)表達(dá)式 x0是根的初值 輸出參數(shù)中 x是非線性方程的數(shù)值解 fx是數(shù)值解的函數(shù)值 返回參數(shù)flag 0時 表示求解成功 否則求解出現(xiàn)問題 函數(shù)fzero所使用的算法為二分法 secant法和逆二次插值法的組合 例8 4求解一維非線性方程 求解單變量x非線性方程x0 0 1 解的初值 xz fz flag fzero atan x exp x x0 disp 求解成功性判斷參數(shù) flagdisp 非線性方程的解 xzdisp 非線性方程解的函數(shù)值 fzM文件運(yùn)行結(jié)果 求解成功性判斷參數(shù) flag 1非線性方程的解 xz 0 6066非線性方程解的函數(shù)值 fz 1 1102e 016 2 多維非線性方程組對于多維非線性方程組使用多維非線性方程組優(yōu)化解函數(shù)求解的格式 x fval flag fsolve fun x0 其中 輸入?yún)?shù)中fun是非線性方程組的向量函數(shù)表達(dá)式 x0是根的初值 輸出參數(shù)中x是非線性方程 組 的數(shù)值解 fval是數(shù)值解的函數(shù)值 返回參數(shù)flag 0時 表示求解成功 函數(shù)fsolve的作用是從根的初值x0開始 以逐漸減少誤差的算法 搜索出滿足多維非線性方程組fun的實(shí)根x和對應(yīng)的函數(shù)值fval 例8 5求解二維非線性方程組x0 2 3 解的初值 定義非線性方程組表達(dá)式f和向量xfun inline 2 x 1 x 2 2 exp x 1 x 1 3 x 1 x 2 exp x 2 x xn fval flag fsolve fun x0 disp 求解成功性判斷參數(shù) flagifflag 0disp 方程組的解成功 elseifflag 0disp 方程組的解不成功 enddisp 非線性方程組的解 xndisp 非線性方程組解的函數(shù)值 fval M文件運(yùn)行結(jié)果 求解成功性判斷參數(shù) flag 1方程組的解成功非線性方程組的解 xn 0 99781 2755非線性方程組解的函數(shù)值 fval 1 0e 006 0 1945 0 3372兩個非線性方程解的函數(shù)值非常接近于0 說明其解的誤差很小 且判斷參數(shù)flag 0 解成功 8 2常微分方程求解 求解微分方程必須事先對自變量的某些值規(guī)定出函數(shù)或是導(dǎo)數(shù)的值 若在自變量為零的點(diǎn)上 給出初始條件 稱為初值問題 最普遍的自變量是 時間 例如 彈性系統(tǒng)的自由振動 若以時間為零來限定位移和速度 這是一個初值問題 若在自變量為非零的點(diǎn)上 給出邊界條件 稱為邊值問題 最普遍的自變量是 位移 例如 描述梁彎曲變形的微分方程 邊界條件總是規(guī)定在梁的兩端 8 2 1常微分方程的解析解在MATLAB中 用大寫的字母D表示導(dǎo)數(shù) Dy表示一階導(dǎo)數(shù) D2y表示二階導(dǎo) 表示微分方程 在MATLAB中 由函數(shù)dsolve進(jìn)行常微分方程 組 的求解問題 其調(diào)用格式是 y dsolve e1 e2 en c1 c2 cn v1 v2 vn 它用于求解常微分方程組e1 e2 en在初值條件c1 c2 cn下的特解 如果不給出初值條件 則求出解常微分方程組的通解 v1 v2 vn是求解變量 例9 1 求常微分方程的通解在此問題中 微分方程的MATLAB表達(dá)式為 Dy y 1 y 2 0沒有給出初值條件 只能求出解常微分方程組的通解 調(diào)用求解微分方程組的語句為 y dsolve Dy y 1 y 2 0 x 運(yùn)算結(jié)果 y 1 1 exp 2 x C1 1 2 1 1 exp 2 x C1 1 2 即該微分方程的解析解為 求常微分方程組的通解 沒有給出初值條件 只能求出解常微分方程組的通解 調(diào)用求解微分方程組的語句為 x y dsolve Dx y exp z Dy exp 2 z z 運(yùn)算結(jié)果 x cos z C2 sin z C1 1 5 exp 2 z 1 2 exp z y sin z C2 cos z C1 2 5 exp 2 z 1 2 exp z 其中 C1與C2是積分常數(shù) 即該微分方程的解析解為 求常微分方程組當(dāng)和條件下的特解 在此問題中 兩個微分方程的MATLAB表達(dá)式為 e1 Dx 2 x Dy 10 cos t e2 Dx Dy 2 y 4 exp 2 t 初值條件表達(dá)式為 C1 x 0 2C2 y 0 0 系統(tǒng)函數(shù)dsolve默認(rèn)自變量為t 所以調(diào)用求解微分方程組的語句為 x y dsolve Dx 2 x Dy 10 cos t Dx Dy 2 y 4 exp 2 t x 0 2 y 0 0 運(yùn)算結(jié)果 x 2 exp t sin t 4 cos t 3 sin t 2 exp 2 t y 2 exp t cos t sin t 2 cos t 即該微分方程的解析解為 8 2 2常微分方程的數(shù)值解MATLAB中用于求解一階微分方程組初值問題數(shù)值解的常用函數(shù)是ode45 它是基于四階五級Runge Kutta變步長算法 采用誤差作為監(jiān)測指標(biāo)的閉環(huán)控制 既保證有較高的運(yùn)算速度和中等計(jì)算精度 也保證預(yù)期的數(shù)值穩(wěn)定性 是大部分場合求解非剛性常微分方程數(shù)值解的首選方法 求解非剛性常微分方程數(shù)值解的常用函數(shù)還有ode23 它是基于二階三級Runge Kutta的單步算法 計(jì)算精度較低 函數(shù)ode45的調(diào)用格式是 t x ode45 fun t0 tf x0 其中 fun是描述微分方程的符號函數(shù)或是用inline 函數(shù)描述的微分方程 t0表示起始時刻 tf表示終止時刻 x0是初值問題的初始狀態(tài)變量 例9 2試求洛倫茨 lorenz 系統(tǒng)動態(tài)模型狀態(tài)方程的數(shù)值解 初值 終止時刻 其中 系數(shù) 描述洛倫茨 lorenz 系統(tǒng)動態(tài)模型狀態(tài)方程3維非線性常微分方程組的函數(shù)文件functionxdot lorenzeq t x flag beta rho sigma 輸入?yún)?shù)中 beta rho sigma是方程組有關(guān)項(xiàng)的系數(shù) 變量flag用于占位xdot beta x 1 x 2 x 3 rho x 2 rho x 3 x 1 x 2 sigma x 2 x 3 求解洛倫茨 lorenz 系統(tǒng)3維非線性常微分方程組tf 100 x0 0 0 1e 10 tf是終止時刻 beta 8 3 rho 10 sigma 28 給各系數(shù)賦值 t x ode45 lorenzeq 0 tf x0 beta rho sigma subplot 1 2 1 plot t x 繪制2維線圖gridontitle bf狀態(tài)變量的時間響應(yīng)圖 subplot 1 2 2 plot3 x 1 x 2 x 3 繪制3維線圖axis 1042 2020 2025 title bf狀態(tài)變量的相空間三維圖 8 2 3高階微分方程的降階求解n階微分方程一般形式 已知輸出變量各階導(dǎo)數(shù)初值選擇一組狀態(tài)變量可以把n階微分方程化成等價(jià)的一階微分方程組這是一個含n個未知函數(shù)的一階微分方程組 初值 例9 3試求二階微分方程的數(shù)值解 狀態(tài)變量初值和 起始時刻 終止時刻 該二階微分方程可以通過下面的變換 得到等價(jià)的一階方程組 編制如下函數(shù)m文件 functiondy weifen t x dy zeros 2 1 dy 1 x 2 x 2 zdy 2 sin t x 2 x 1 x 1 y 然后用函數(shù)ode45求解微分方程 t y ode45 weifen 0 20 0 6 plot t y 1 holdon 繪制y的圖像plot t y 2 r 繪制y 的圖像xlabel bft ylabel bfy title bf二階微分方程的數(shù)值解 legend 實(shí)線 y的圖像 虛線 Dy的圖像 運(yùn)算結(jié)果如圖8 3所示 例9 4試求二階微分方程的數(shù)值解 該二階微分方程可以通過下面的變換 得到等價(jià)的一階方程組 已知狀態(tài)變量的初值起始時刻 終止時刻 1 計(jì)算二階微分方程的解析解symsty dsolve D2y 3 cos 2 t 2 sin t t 3 5 y 0 0 Dy 0 0 t ezplot y 010 繪制解析解xlabel bft ylabel bfy title bf二階微分方程的解析解和數(shù)值解 gtext bfy 3cos2t 4 2sint t 3 6 7t 2 4 2t 3 4 holdon 2 計(jì)算二階微分方程的數(shù)值解fun inline x 2 3 cos 2 t 2 sin t t 3 5 t x tx1 ode45 fun 0 10 0 0 plot t x1 1 r 繪制數(shù)值解gridonlegend 曲線 解析解 星號 數(shù)值解 運(yùn)算結(jié)果如圖8 4所示 8 2 4剛性微分方程的數(shù)值解剛性方程是指其Jacobian矩陣的特征值相差十分懸殊 在解的性態(tài)上表現(xiàn)為 其中一些解變化緩慢 另一些變化快 且相差較懸殊 這類方程常常稱為剛性方程 又稱為Stiff方程 剛性方程和非剛性方程對解法中步長選擇的要求不同 剛性方程一般不適合由ode45這類函數(shù)求解 而應(yīng)該采用ode15s等函數(shù) 如果不能分辨是否是剛性方程 先試用函數(shù)ode45 再用函數(shù)ode15s 求解剛性微分方程數(shù)值解常用函數(shù)有 ode15s 它是基于Gears反向數(shù)值微分的多步算法 具有中等計(jì)算精度 ode23s 它是基于Rosebrock的單步算法 計(jì)算精度較低 函數(shù)ode15s和ode23s的調(diào)用格式與函數(shù)ode45相同 例9 5求剛性微分方程組的數(shù)值解已知 初始條件和 起始時刻和終止時刻 使用函數(shù)dsolve可求出解析解為 可見 當(dāng)時 兩個分量均趨于0 下降極快 而下降很慢 編制描述剛性微分方程組的M文件 剛性微分方程組的函數(shù)文件functionfun stiff t y fun 0 01 y 1 99 99 y 2 100 y 2 調(diào)用函數(shù)ode15s求解 t y ode15s stiff 0 400 2 1 tstep length t 總步數(shù)minh min diff t 最小步長maxh max diff t 最大步長plot t y xlabel bft ylabel bfy title bf剛性微分方程組的數(shù)值解 gtext bfy 1 gtext bfy 2 8 2 5微分代數(shù)方程的數(shù)值解微分代數(shù)方程 differentialalgebraionequation DAE 是指微分方程中某些變量受到某些代數(shù)方程條件的約束 它的一般形式是 其中 描述的方法與普通的微分方程相同 為奇異矩陣 由M文件的求解選項(xiàng)變量Mass說明 在無法嚴(yán)格證明是奇異矩陣時 可以先假設(shè)它是奇異矩陣 如果在MATLAB求解過程中沒有出現(xiàn)是非奇異矩陣的警告信息 則求解是有效的 否則無效 例9 6求微分代數(shù)方程組的數(shù)值解 已知 初始條件和 起始時刻和終止時刻 該方程組的第3個是代數(shù)方程 給定3個狀態(tài)變量之間的約束關(guān)系 用矩陣形式表示該微分代數(shù)方程 編寫描述微分代數(shù)方程的函數(shù)文件 functiondx DAE t x dx 0 2 x 1 x 2 x 3 0 3 x 1 x 2 2 x 1 x 2 5 x 2 x 3 2 x 2 x 2 x 1 x 2 x 3 1 Options是用命令odeset設(shè)置的可選積分參數(shù) 將常系數(shù)矩陣M賦給質(zhì)量矩陣Mass M 1 0 0 0 1 0 0 0 0 x0 0 8 0 1 0 1 options odeset options Mass M t x ode15s DAE 0 20 x0 options plot t x 1 r t x 2 b t x 3 g title 微分代數(shù)方程的數(shù)值解 gtext bfx 1 t gtext bfx 2 t gtext bfx 3 t xlabel bft ylabel bfx gridon